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ABSTRACT 

This  thesis  deals  with  two  closely  related  aspects 
of  the  study  of  non-homogeneous  seepage.  Part  I  undertakes 
the  calculation  of  pressure  distributions  using  the  finite 
element  method  and  Part  II  describes  attempts  to  measure 
water  pressures  with  a  new  piezometer  system. 

In  Part  I  a  technique  is  developed  and  utilized 
vzhich  takes  into  account  the  effects  of  stress  dependant 
permeability  in  a  rock  material  upon  the  pressure  distri¬ 
bution  in  unconfined  flow  situations.  Finite  element  seepage 
theory  is  reviewed  and  the  finite  element  seepage  method  is 
coupled  with  the  finite  element  elastic  stress  analysis  to 
perform  the  subsequent  computations.  A  study  of  permea¬ 
bility  measurements  has  been  undertaken  and  a  function 
relating  permeability  to  stress  has  been  derived.  This 
function  is  used  in  the  analyses. 

A  configuration  of  a  steep-walled  excavation  is 
studied  with  variations  in  initial  stress  conditions  and 
permeability  functions.  Equipotential  plots,  selected 
contours  of  isotropic  permeability  and  contours  of  average 
principal  effective  stress  are  presented  and  discussed. 

The  factors  that  affect  the  pressure  distributions  are  noted. 

Part  II  deals  with  the  field  installation  of  the 
University  of  Alberta  multistage  piezometer  system  in  the 
Upper  Cretaceous  Edmonton  Formation.  Installation  and 
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measurement  procedures  are  discussed.  Perched  water 
tables  are  shown  to  exist.  Recommendations  and  conclu¬ 
sions  are  formed. 
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CHAPTER  I 


INTRODUCTION 


1. 1  General 

The  solution  of  non-homogeneous  seepage  problems 
is  of  concern  in  many  circumstances.  The  prediction  of 
flow  quantities  or  pressure  magnitudes  is  an  important 
consideration  in  geotechnical,  water  resource  and  irriga¬ 
tion  engineering.  Due  to  the  many  combinations  and  varia¬ 
tions  that  occur  in  the  permeabilities  of  non-homogeneous 
media,  a  great  diversity  of  flow  regimes  exist. 

In  soil  and  rock  mechanics  the  knowledge  of  induced 
or  natural  groundwater  conditions  is  of  paramount  conse¬ 
quence  in  stability  analyses.  A  familiarity  with  the 
factors  which  influence  these  conditions  is  then  necessary 
in  order  to  properly  evaluate  the  most  critical  conditions 
for  instability  to  occur.  The  variation  of  permeability 
throughout  the  media  is  an  important  factor. 

1. 2  The  Nature  of  Non-homogeneous  Flow  in  the  Field 

In  order  to  fully  appreciate  the  complexities  of 
non-homogeneous  seepage  problems,  the  overall  setting  in 
which  flow  takes  place  has  to  be  examined.  Groundwater 
movement  can  be  viewed  initially  in  the  framework  of  the 
hydrological  cycle  that  occurs  in  nature.  It  may  then  be 
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seen  in  the  context  of  the  topographical  and  geological 
formation  of  the  area.  With  this  perspective  the  ground- 
water  conditions  in  a  valley  or  excavation  with  their  own 
unique  stratigraphy  and  geometric  configurations  are  more 
fully  understood.  Further  complications  will  occur  when 
the  seepage  conditions  caused  by  draining  or  damming  are 
superimposed  on  the  "natural"  flow  regime. 

The  permeabilities  of  the  materials  in  a  geological 
formation  and  also  on  the  smaller  scale  of  an  individual 
slope  may  be  quite  varied  in  terms  of  both  non-homogeneity 
and  anisotropy.  This  in  turn  will  affect  the  groundwater 
flow  characteristics  of  the  area.  Impermeable  strata 
situated  among  more  permeable  layers  often  result  in  perched 
water  tables.  For  example  thin  layers  of  bentonite  of  low 
permeability  that  are  underlain  by  a  draining  coal  seam  tend  to 
produce  this  condition.  In  other  areas  porous  limestone 
resting  on  impermeable  shale  may  bring  about  the  same 
result.  A  fault  in  a  rock  slope  may  act  as  a  groundwater 
barrier,  groundwater  conduit  or  as  a  subsurface  drain.  A 
slope  may  be  comprised  of  a  relatively  homogenous  rock 
material  but  irregular  jointing  may  cause  large  variations 
in  the  equivalent  permeabilities  of  the  joints  resulting  in 
irregular  pressure  distributions  in  the  slope.  These 
examples  coupled  with  a  realization  of  the  effects  of 
regional  groundwater  behaviour  (Toth,  1962;  Freeze  and 
Wither spoon ,  1966,  1967,  1968)  are  sufficient  to  illustrate 
the  vast  differences  that  can  occur  in  nature  in  groundwater 
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Various  workers  have  reported  cases  and  situations 
where  the  consideration  of  non-homogeneous  flow  has  been 
essential.  Lyell  (1969)  mathematically  formulates  the 
problem  where  several  layers  drain  into  an  open  pit  mine 
resulting  from  a  number  of  perched  water  tables.  He  uses 
field  piezometer  measurements  in  conjunction  with  the 
analytical  studies  to  predict  future  seepage  patterns.  With 
this  information  the  effects  of  slope  stability  are  evaluated. 
In  another  case,  Steffen  and  Klingman  (1966)  indicate  that 
piezometric  data  is  invaluable  for  monitoring  and  evaluating 
the  stability  of  slopes  at  the  Nchanga  open  pit  mine. 

While  arriving  at  shear  strength  values  in  a  failed  slope, 
Eigenbrod  and  Morgenstern  (1971)  illustrate  the  importance 
of  interpreting  the  water  pressure  distribution  correctly 
at  a  location  where  several  perched  water  tables  occur. 

Khan  (1970)  attempts  to  show  that  the  effects  of  non- 
uniform  permeabilities  aid  in  explaining  in  part  the  Vajont 
reservoir  disaster.  In  a  type  of  retrogressing  slope 
failure  problem,  Henkel  (1967)  has  described  the  effect  of 
flow  in  non-homogeneous  deposits  and  shoves  how  under 
certain  conditions  degradation  of  the  slope  occurs  due  to 
its  non-homogeneity. 

Generally  the  effects  of  foundation  permeability 
variations  are  an  essential  consideration  in  dam  design. 

One  case  where  these  variations  were  assumed  to  be  relatively 
inconsequential  illustrates  dramatically  the  importance  of 
non-homogeneous  flow.  In  the  Malpasset  dam  failure  (Bernaix, 
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1969) ,  the  effect  of  confining  stresses  imposed  by  the 
structure  on  the  rock  foundation  material  caused  greater 
permeability  variations  than  anticipated  with  the  result 
that  unexpected  critical  seepage  gradients  led  to  the 
catastrophe . 

These  examples  should  suffice  to  show  the  importance 
of  predicting  as  well  as  measuring  the  seepage  behaviour  in 
non-homogeneous  media.  For  a  more  complete  and  integrated 
understanding  of  non-homogeneous  flow  it  is  evident  that 
both  analysis  and  instrumentation  will  continue  to  play  an 
important  role  in  the  future. 

1 . 3  Statement  of  the  Problems 

Recently  seepage  and  groundwater  studies  have 
received  renewed  interest  since  the  digital  computer  has 
facilitated  problem  solving.  Until  recently  seepage 
analyses  of  non-homogeneous  media  have  been  difficult  to 
handle.  Various  attempts  using  theoretical,  experimental 
and  analogue  methods  have  usually  resulted  in  cumbersome 
and  inflexible  means  of  solution. 

Computer  applications  using  finite  element  and 
finite  difference  numerical  procedures  may  be  used  in  the 
treatment  of  non-homogeneous  seepage  problems  (Zienkiewicz 
and  Cheung,  1965;  Tomlin,  1965).  Finite  element  methods 
have  proved  to  be  the  most  flexible  of  these  numerical 
analyses  and  it  is  with  this  method  that  unconfined  flow 
problems  have  also  been  solved.  Early  finite  element 


seepage  programs  did  not  include  the  solution  of  these 
problems  due  to  the  boundary  conditions  that  are  present 
at  the  free  surface,  i.e.,  the  velocity  normal  to  the  free 
surface  vR  =  0,  and  the  pressure  at  the  free  surface,  P  =  0. 
Various  workers  (Taylor  and  Brown,  1967;  Finn,  1967;  Neuman 
and  Witherspoon,  1970)  have  incorporated  these  boundary 
conditions  in  conventional  seepage  programs  through  iteration 
techniques  using  mesh  and  nodal  movement. 

With  the  finite  element  method  it  is  also  possible 
to  solve  steady  state  free  surface  flow  in  a  medium  of 
continuously  varying  isotropic  or  anisotropic  permeabilities. 
One  important  case  of  the  latter  comes  about  when  the 
permeabilities  are  dependent  on  the  imposed  state  of  stress. 
This  may  occur  in  rock  slopes  where  the  secondary  permeability 
(i.e.  along  the  cracks  or  joints)  is  the  major  source  of 
water  conductivity.  This  phenomenon  of  permeability  changes 
with  stress  has  been  observed  in  the  laboratory  (e.g. 

Patching,  1965;  Bernaix,  1969)  and  a  permeability  decrease 
with  depth  which  is  attributed  in  part  to  increased  stress 
has  been  observed  in  the  field  (e.g.  Davis  and  Turk,  1964; 
Sherman  and  Banks,  1970).  For  example,  from  the  work  of 
Bernaix  and  Patching  it  is  known  that  stress  changes  in 
samples  of  various  materials  such  as  fissured  gneisses  and 
friable  coals  cause  changes  in  permeability  that  may  be 
several  orders  of  magnitude  depending  on  the  extent  of  the 
stress  changes.  One  of  the  problems  of  concern  here  was 
to  analyze  the  resultant  pressure  distribution  in  an 
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excavation  where  the  permeability  throughout  the  medium  was 
a  function  of  the  stress. 

When  determining  the  flow  behaviour  in  a  slope,  it 
is  obvious  that  as  the  degree  of  non-homogeneity  increases, 
the  need  for  adequate  piezometric  instrumentation  also 
increases  (e.g.  Steffan  and  Klingman,  1966) .  However  it 
is  desirable  to  make  an  estimate  of  the  groundwater  condi¬ 
tions  in  the  most  economical  and  efficient  manner.  Often 
a  multitude  of  measuring  elevations  and  locations  are 
required  (e.g.  Lyell,  1969)  and  to  facilitate  measurements, 
it  is  convenient  to  utilize  a  piezometer  with  rapid  response 
characteristics  (Penman,  1961).  Electric  diaphragm  piezo¬ 
meters  incorporate  such  a  feature  (Brooker  et  al,  1968) 
and  the  pressure  measurement  itself  can  be  determined  with 
great  ease.  Unfortunately  drawbacks  of  diaphragm  piezometers 
in  general  have  been  their  lack  of  dependability  (Lindberg, 
1965)  as  well  as  the  excessive  costs  of  electrical  diaphragm 
tips.  Hence  in  addition  to  taking  into  account  their 
dependability,  it  is  desirable  to  decrease  these  costs.  As 
a  result  piezometer  tips  have  been  developed  which  are 
recoverable  and  initially  these  instruments  were  designed 
to  measure  the  pressure  at  one  location  per  borehole 
(Lundgren,  1966).  In  a  second  task  this  thesis  assesses 
the  University  of  Alberta  multistage  piezometer  system. 

This  electrical  diaphragm  piezometer  is  designed  to  measure 
pressures  at  various  pre-selected  levels  in  a  borehole. 

With  the  pressure  information  obtained  one  might  also 
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evaluate  the  effects  of  non-homogeneity  on  the  seepage 
behaviour . 


PART  ONE 


CHAPTER  II 

ANALYSIS 


2 . 1  Introduction 

This  chapter  discusses  the  finite  element  method 
of  analysis  of  problems  in  steady  seepage  and  related 
aspects  of  concern.  A  seepage  program  (Taylor  and  Brown, 
1967)  is  coupled  with  a  finite  element  linear  elastic 
stress  program  to  solve  for  the  phreatic  surface  and 
pressure  distribution  in  a  medium  with  stress  dependent 
isotropic  permeability.  Taylor  and  Brown's  program  is 
used  to  find  the  seepage  body  forces  for  the  stress  analysis. 
After  the  body  forces,  calculated  from  the  seepage  distri¬ 
bution,  have  been  assigned  throughout  the  stress  grid  the 
effective  stresses  are  computed  and  used  to  find  the 
permeabilities.  These  are  taken  to  be  a  function  of  the 
effective  stresses  in  the  manner  described  in  Chapter  III. 

The  permeabilities  are  assigned  to  the  proper  seepage 
elements  and  the  stable  position  of  the  free  surface  is 
found  by  iteration.  This  cycle  of  finding  the  body  forces, 
stress  distribution  and  tree  surface  location  is  repeated 
until  the  difference  in  the  stable  positions  of  the  free 
surface  of  two  successive  cycles  indicates  that  convergence 
is  obtained.  The  analysis  is  applied  to  seepage  into  an 

N 

excavation  as  illustrated  in  Chapter  IV. 


8 


2.2 
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Derivation  of  the  Seepage  Equation 

The  flow  of  water  in  a  porous  medium  is  caused  by 
a  hydraulic  gradient,  i,  which  may  be  defined  as 


i  =  gradient  (Y  +  P/yw)  =  grad  U 


(2.1) 


where  Y  =  level  at  the  point  related  to  some  datum 
P  =  pressure  at  the  point 
yw=  unit  weight  of  water 
U  =  hydraulic  potential 

Darcy's  law  states  that  the  apparent  velocity,  v,  is 
proportional  to  the  gradient  i  of  the  hydraulic  potential. 


v  =  -ki  =  -k  grad  U 


(2.2) 


If  the  anisotropic  nature  of  the  medium  is  considered 
in  two  dimensions  the  constant  of  proportionality  known  as 
the  permeability,  k,  is  given  by  a  matrix  of  values  and  may 
be  written  as, 


^xx  kxy 


kyx  kyy 


(2.3) 


and  therefore  the  velocity  v  is  a  vector  given  by 
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v  1 
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<  x  =  - 

XX 

xy 
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V  1 

k 

k 

%  U 
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_  yx 
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l3y  i 

(2.4) 


If  x  and  y  are  the  principal  directions  of  aniso¬ 
tropy  and  k  and  k  are  the  principal  values  of  permeability, 
x  y 

the  velocities  may  be  expressed  as 


k 

x  3x 


(2.5a) 


v 


v 


(2.5b) 


Imposing  the  condition  of  continuity  upon  a  unit 
volume  and  assuming  an  incompressible  fluid,  the  following 
is  obtained. 


9  vv  9  v.. 
— x  +  — I 

8  x  9  y 


0 


From  (2.5)  this  becomes. 


9  /  9U  9  /.  9  U  \ 

7“  i  kx —  ,  +  — (  kv —  =  0 

9x  '  x9 x  /  9  y  \  y9  y  / 


(2.6) 


(2.7) 


In  this  work. 


k  = 


=  k  (x ,  y ) 


(2.8) 
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and  hence  equation  (2.7)  becomes. 


9k  3  U  9  k  9U  k3  2U  k  92U 

- + - +  — y  ~  =  0 

9x  8  x  9  y  ay  9xz  d  yz 


2 . 3  Finite  Element  Formulation 

Even  though  isotropic  permeabilities  are  used  in 
this  work,  anisotropy  will  be  considered  in  the  following 
formulation  of  the  finite  element  analysis. 

A  permeability  matrix  in  the  x-y  global  system' is 
defined  in  terms  of  the  permeabilities  which  are  usually 
determined  parallel  and  perpendicular  to  the  principal 
axis  of  permeability.  The  axes  may  be  at  an  angle  0  to 
the  global  x  axis  as  shown  in  Figure  2.1. 

Quantities  in  the  local  co-ordinate  system  are 

primed . 

From  equations  (2.5),  velocities  in  the  local 
system  may  be  expressed  as 


/ 
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0 
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(2.10) 


It  can  easily  be  shown  that 
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(2.11) 
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FIG.  2.1.  LOCAL  AND  GLOBAL  CO-ORDINATE  SYSTEMS 
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and 

(  3^0 
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From  these  equations  and  equation  (2.4)  the  permea¬ 
bility  matrix  in  the  global  system  may  be  related  to  the 
values  in  the  local  system  as  follows. 


k 


(2.13) 


To  use  the  finite  element  method  the  properties  of 
an  individual  element  have  firstly  to  be  established.  The 
following  derivation  will  be  for  the  element  known  as  the 
constant  velocity  triangle. 

The  simplest  assumption  for  the  distribution  of 
pressures  P  across  the  element  in  Fig.  2.2a  is  that  they 
vary  linearly  between  the  vertices  (Zienkiewicz  and  Cheung, 
1965;  Zienkiewicz,  Mayer  and  Cheung,  1966).  Hence  the 
pressure  P  at  any  location  throughout  the  triangle  will  be 
given  by 


or 


P  =  +  c*2x  +  «3y 


a2  >■ 

Lj 


P  = 


X 


(2.14a) 


(2.14b) 


The  three  unknowns  in  equation  (2.14)  may  be 
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FIG.  2.2b.  ELEMENT  FLOWS 


FIG.  2.2 


THE  SEEPAGE  ELEMENT 
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evaluated  by  forming  equations  in  o^,  a0  and  a3  at  the 
three  vertices 


f  D  _  " 


={  P2  = 

1  ^  i 


^P3' 


1  xi  yi 

(Vl 

1  x2  y2 

a2  > 
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Therefore 
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and 


p  =  [i x  yj  [A_1 
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where 

r— 

tx  t2  t3 

E*'1] 

1 

y23  y31  y12 

tl  +  t2  +  t3 

Lx32  x13  x 2 1 J 

(2 
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and  x . .  =  x . -x . 
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Vij  =  yi-yj 


Pressure  gradients  are  obtained  as 
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Using  equations  (2.4)  and  (2.19)  velocities  may  be 


given  as 
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where 


P  + 


/ill 
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(2.20) 


BY 
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V 

and  P  will  have  units  of  length 


Equivalent  nodal  point  flows  may  be  written  for 
each  element.  Then  continuity  of  flow  at  node  N  requires 
that 


£ 

E  V  =  0 

i=l  (2.21) 


where  %  is  the  number  of  elements  connected  to 
node  N 


From  Fig.  2.2b  the  horizontal 
flow  Qy  for  node  i  may  be  written  as 


flow  Qx  and  vertical 


(2.22a) 


(2.22b) 


. 
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by 


Hence  the  nodal  flows  may  be  related  to  the  velocities 


(q)  =  I  L 


'  ^ 
jV'i 
> 


(2.23) 


where 


= 


23 

x32 

31 

x13 

12 

X21 

(2.24) 


Using  equation  (2.20),  (2.23)  becomes 
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+  L  !!  k 


1 1 

JL 
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I8Y 

a  S 


(2.25) 


Equation  (2.25)  is  used  in  forming  equation  (2.21). 
Computing  equation  (2.21)  requires  solving  a  set  of  simul¬ 
taneous  equations  in  terms  of  P  .  The  nodes  where  the 
pressures  are  specified  will  have  to  be  modified  to  give 
the  known  pressure  as  the  solution.  Equation  (2.9)  is  the 
equation  of  steady  state  flow  which  is  solved  for  the  two 
dimensional  case  assuming  certain  boundary  conditions  by 
equation  (2.21)  through  the  assignment  of  the  appropriate 
permeability  to  each  discrete  element.  Equation  (2.21) 
has  been  derived  considering  an  assembly  of  elements  while 
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equation  (2.9)  was  obtained  assuming  a  continuum.  Both 
have  been  deduced  by  giving  heed  to  continuity  and  Darcy's 
law. 

The  program  of  Taylor  and  Brown  can  use  either 
triangles  or  quadrilaterals.  When  quadrilateral  elements 
are  used  they  are  divided  into  four  triangular  elements  by 
two  diagonals.  A  provision  is  included  for  generation  of 
quadrilateral  elements.  Additional  details  for  the  use  of 
the  program  are  given  by  Kealy  and  Busch  (1970) . 

A  similar  derivation  corresponding  to  the  one 
presented  above  may  be  given  for  the  stress  analysis. 
Instead  of  pressures  at  the  nodes  displacements  are  used 
in  the  displacement  formulation.  Since  displacement  is 
a  vector  quantity  with  two  components  for  a  two  dimensional 
analysis, the  matrices  are  expanded  but  the  same  procedure 
may  be  followed  as  outlined  here.  Holand  and  Beil  (1970) 
give  details  of  the  method. 

2 . 4  Free  Surface  Location  Determination 

The  free  surface  is  located  where  the  flow,  vR 
normal  to  the  surface  is  zero  and  the  pressure,  P,  at  the 
surface  is  atmospheric.  Due  to  the  number  of  equations 
involved  either  vn  can  be  set  to  zero  or  P  can  be  set  to 
zero.  The  program  of  Taylor  and  Brown  sets  vn  =  0  and 
iterates  so  that  P  -*»  0 . 

During  an  iteration  the  solution  of  the  finite 
element  equations  yields  positive  or  negative  pressures. 


p  u 

■ 


19 


P,  at  the  free  surface  nodes.  In  the  program  a  location 
correction  for  any  free  surface  node  M  is  made  according 
to  the  equation 

dx  =  BETA  x  p  x  cos0/yw  (2.26a) 

dy  =  BETA  x  p  x  sin0/yw  (2.26b) 

where  dx  and  dy  are  the  correction  components, 

BETA  is  a  factor  to  prevent  instability  of 
the  iteration  procedure,  and 

0  is  the  correction  direction  with  respect 
to  the  x  axis  for  node  M. 

If  the  nodes  are  corrected  to  where  the  iteration 
gives  zero  excess  pressure,  then  the  position  may  be 
overcorrected  (i.e.  if  BETA  =  1.0)  and  instability  may 
occur.  To  prevent  this  a  value  of  BETA  <1  is  assigned. 
Generally  a  value  of  0.85  to  0.90  is  used.  Fig.  2.3 
demonstrates  the  procedure. 

The  program  calculates  the  square  root  of  the  sum 
of  the  squares  of  the  x  and  y  correction  components  for  n 
free  surface  nodes  and  compares  this  to  the  input  tolerance 
value  TOL.  If 

TOL  >  J dxi2  +  dyi2  +  dx22  +  dy22  +  •••  +  dxn2  +  dYn2 

(2.27) 

the  iteration  process  is  stopped.  A. tolerance  of  0.1%  of 
the  net  head  may  be  used.  In  general  if  the  free  surface 
is  converging  properly  a  tolerance  of  0.01%  will  produce 
negligible  difference  in  the  final  result  compared  to  a 
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FIG.  2.3 


FREE  SURFACE  NODE  CORRECTION  PROCEDURE 
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tolerance  in  the  order  of  0.1%. 

The  program  generates  nodal  points  by  spacing  them 
evenly  on  a  straight  line  between  any  two  consecutive  input 
nodes  with  a  nodal  number  difference  greater  than  one. 

After  adjustment  of  the  free  surface  nodes  the  intermediate 
points  are  regenerated  and  hence  the  nodes  between  the  free 
surface  and  the  adjacent  input  nodes  will  be  "strung"  along. 

In  order  to  gain  confidence  in  the  procedure  it  is 
best  that  the  discretization  of  the  finite  element  method 
be  checked  with  a  closed  form  free  surface  solution.  A 
closed  form  solution  as  computed  by  Gunther  (Mu skat ,  1946) 
for  the  phreatic  line  of  a  configuration  that  might  occur 
in  the  core  of  a  rockfill  dam  is  shown  as  a  dotted  line 
superimposed  on  the  finite  element  solution  in  Fig.  2.4. 

It  may  be  seen  that  there  is  a  close  correspondence  between 
the  two  solutions  although  the  finite  element  solution  is 
slightly  lower  for  most  of  the  free  surface  and  slightly 
above  it  at  the  exit  face.  When  the  free  surface  location 
of  a  30°  slope  was  checked  with  the  finite  difference  solu¬ 
tion  of  Khan  (1971)  the  same  sort  of  behaviour  occurred  but 
the  discrepancy  for  some  parts  of  the  free  surface  appeared 
slightly  larger  in  terms  of  percentage  of  free  surface 
height  from  the  base.  However,  Khan's  comparison  with  the 
free  surface  solution  of  Numerov  (Harr,  1962)  shows  the 
finite  difference  solution  to  be  somewhat  higher.  A  close 
approximation  to  the  solution  at  the  face  in  Fig.  2.4  is 
difficult  since  the  free  surface  approaches  it  tangentially 
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FINITE  ELEMENT 

COMPARISON  WITH  CLOSED  FORM  SOLUTION 
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- GUNTHER'S  CLOSED  FORT!  SOLUTION 


FIG.  2.4 
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(Casagrande ,  1940)  .  Also,  Taylor1  and  Brown  (1967)  stats , 

An  ambiguity  .  .  .  must  exist  where  the  free  surface  inter— 

sects  the  geometric  boundary  of  the  body  because  the  nodal 
flow  at  this  point  cannot  be  zero.  This  difficulty  is 
minimized  by  reducing  the  mesh  size  adjacent  to  this  point." 

In  some  cases  a  closer  approximation  may  be  obtained  by 
taking  a  greater  frequency  of  elements  near  the  face. 

However  this  may  create  a  second  problem  if  in  the  final 
grid  shape  an  abundance  of  thin  elongate  odd  shaped  triangular 
or  quadrilateral  elements  occur.  This  will  usually  result  in 
ill- conditioning  of  the  finite  element  equations  causing  a 
lack  of  convergence  or  giving  the  improper  location  of  the 
free  surface.  The  second  problem  may  be  rectified  by  esti¬ 
mating  the  final  position  of  the  free  surface  and  designing 
a  grid  with  both  triangular  and  quadrilateral  elements  if 
required  so  that  the  final  grid  configuration  will  contain 
a  minimum  of  odd  shaped  elements. 

Ill-conditioning  occurs  when  small  changes  in  the 
coefficients  of  a  linear  system  cause  large  changes  in  the 
solution.  Small  errors  will  occur  in  the  banded  flow7  matrix 
during  its  solution  because  of  the  accumulation  of  round¬ 
off  error.  It  was  discovered  that  ill-conditioning  occurred 
in  these  problems  when  certain  triangular  or  quadrilateral 
element  shapes  caused  one  term  on  the  main  diagonal  of  the 
banded  flow  matrix  to  be  much  smaller  then  the  rest. 
Unfortunately  an  attempt  at  using  double  precision  computation 
resulted  in  no  apparent  improvement.  Conte  (1965)  discusses 
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the  phenomena  of  ill-conditioning  of  linear  systems  and  may 
be  consulted  for  a  fuller  treatment  of  the  subject. 

The  results  of  poor  element  shape  has  been  noted 
in  the  work  of  Kealy  and  Busch  (1971)  and  Neuman  and  Wither¬ 
spoon  (1970).  This  is  illustrated  in  Figure  2.5a  and  2.5b. 
Neuman  and  Witherspoon  compare  the  results  of  their  more 
rapidly  converging  rinite  element  solution  with  the  results 
of  Taylor  and  Brown's  and  attribute  the  divergence  to  an 
inconsistent  method  of  nodal  shifting  and  fail  to  indicate 
that  a  judicious  choice  of  element  configuration  will  not 
cause  divergence. 

The  author  had  problems  similar  to  these  mentioned 
here.  The  pressure  at  the  second  free  surface  node  from 
the  face  would  not  tend  to  zero  giving  the  sort  of  effect 
shown  in  Fig.  2.5a.  It  was  found  that  the  ratio  of  the 
smallest  to  the  largest  term  of  the  odd  shaped  element  flow 
matrices  was  in  the  order  of  1/200  when  this  occurred. 

When  the  element  shapes  were  modified  in  this  region  as 
shown  in  Fig.  2.6  this  ratio  was  decreased  to  the  order  of 
1/20  and  the  pressure  at  the  second  node  tended  to  zero 
accordingly.  The  author  feels  that  the  effect  of  element 
shape  should  be  investigated  more  fully  on  a  methodical 
basis.  It  might  also  be  advantageous  to  probe  the  possibility 
of  modifying  the  solution  procedure  of  the  banded  flow  matrix 
so  that  ill-conditioning  of  the  finite  element  equations  will 


be  minimized. 
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FFV500  3  SOIL  LOTI  (FINAL)  K1=KP-K3  KH-KV  STRATIFICATION  =  B  DELREfS 
FINITE  ELELENT  LE9t 
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.ELEMENTS  CAUSING  ILL-CONDITIONING 


Output  From  o  Finite-Element  Model  of  a  Homogeneous  Embankment  on  an 
Impermeable  Base. 


FIG.  2.5a  ODD  SHAPED  ELEMENTS  IN  THE  WORK  OF  KEALY  AND 

BUSCH  (1971) 


, —  SOLUTION  AFTER  20 
ITERATIONS  USING 
TAYLOR  AND 
BROWN'S  PROGRAM 


Seepage  through  a  dam  with  toe  drain. 


FIG.  2.5b  ILLUSTRATION  BY  NEUMAN  AND  WITHERSPOON  (1970) 

OF  TAYLOR  AND  BROWN'S  DIVERGING  SOLUTION 


FIG 
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EFFECTS  OF  POOR  ELEMENT  SHAPE 
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FIG.  2.6  MODIFICATION  OF  GRID  CONFIGURATION  TO  PREVENT 

ILL-CONDITIONING  OF  THE  FINITE  ELEMENT  EQUATIONS 
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2.5  Validity  of  the  Results  Using  Varying  Permeability 

For  the  case  of  flow  through  a  horizontal  tube  with 
varying  permeability 

Vx  =  k(x)A£  =  CONSTANT  =  c  (2.28) 


and 


c£x 

Ap  =  (2.29) 

If  an  equal  pressure  drop  across  each  element  is 
desired  Ax/k(x)  must  be  kept  in  the  same  ratio  or 
Ax  a  k(x).  The  smaller  the  permeability,  the  smaller  Ax 
should  be  to  maintain  the  same  ratio.  For  the  two  dimen¬ 
sional  case  and  assuming  k  is  constant, 

v  —  ki  —  kA£  (2.30) 

As 

For  increasing  gradient  and  AP  constant,  As  must  decrease 
or  As  a  Hence  in  general  element  size  should  be 

reduced  in  areas  of  large  gradients  and  low  permeabilities. 

For  the  case  of  flow  through  a  tube  with  permeability 
varying  linearly  in  the  direction  of  flow  a  theoretical  curve 
may  be  obtained  for  the  pressure  along  the  tube  and  compared 
to  the  result  using  various  element  sizes.  However  since  it 
is  difficult  to  evaluate  a  complicated  two  dimensional  case 
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by  extension  of  the  results  from  a  simple  case  it  is  best 
to  observe  the  consequences  with  variation  of  element  size 
in  the  actual  problem. 


2 . 6  Method  of  Analysis 

In  Chapter  III  a  relation  between  the  average 
principal  stress  and  permeability  is  derived.  For  an 
excavation  in  plain  strain  the  average  principal  effective 
stress  may  be  obtained  as  follows: 


The  first  stress  invariant  =  In=  a  +a  +a 

1  x  y  z 


a 


+  a 


I 


II 


III 


(2.31) 


For  the  excavation  process 


AI1  =  AQx  +  Aay  +  A  a  z 


Aa I  +  ^  II  +  A aIII 


(2.32) 


The  final  invariant  of  stress 


iif  =  xi  +  aii 


(2.33) 
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For  plain  strain 

%  =  V(ax  +  ay)  (2.34a) 

and 

A  °z  =  V(A°X  +  &  y)  (2.34b) 

Therefore 

Ilf  =ay(l  +  2K0)  +  (APX  +  A  Oy)  (1  +  v  ) 

(2.35) 

The  average  principal  effective  stress  will  be 

J1  =  k\°y  (1+2K0)  +  (Aax  +  A av)  (1  +v  T 

J  (2.36) 

In  general  for  an  excavation  a  situation  as  shown 
in  Fig.  2.7a  will  be  established.  Since  it  was  thought 
that  the  effects  of  stress  dependent  permeability  would  be 
most  pronounced  for  a  steep  slope  angle,  it  was  decided  to 

study  the  case  of  an  excavation  with  a  vertical  face.  The 

/ 

water  table  was  chosen  to  be  initially  at  the  surface. 

This  simplified  the  problem  since  this  theoretically  pre¬ 
cluded  a  non-triangular  stress  release  distribution  on  the 
face  of  the  excavated  slope. 

Consider  the  following  steps  as  shown  in  Fig.  2.7b- 
2.7e  in  arriving  at  a  method  of  solution.  Fig  2.7b  shows 
the  initial  state  with  no  excavation  and  the  water  table 
at  the  surface.  The  assumed  insitu  stress  conditions  are 
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as  shown.  If  an  excavation  is  made  as  shown  in  Fig.  2.7c 
a  release  of  effective  stresses  will  with  a  finite  element 
program  produce  the  correct  displacements  and  stress  changes 
in  the  slope  provided  that  the  excavation  is  full  of  water. 
If  the  water  is  now  removed  from  the  excavation  a  phreatic 
surface  S  will  develop  as  shown  in  Fig.  2.7d.  Considering 
the  region  above  the  free  surface  marked  II  and  assuming 
the  pressure  to  be  at  atmospheric  above  the  free  surface 
the  change  in  body  force/unit  volume  will  be  equal  to  the 
unit  weight  of  water  since  the  buoyant  force  of  the  water 
is  removed  from  this  area.  In  the  area  marked  I  the 
change  in  body  force/unit  volume  will  be  due  to  the  seepacre 
force  of  the  wTater .  Therefore  to  obtain  the  correct  change 
in  effective  stresses  body  forces  must  be  considered.  Then 
the  final  state  of  effective  stress o ^ '  will  be  given  by 
the  addition  of  the  initial  state,  a  '  the  change  due  to 
stress  release  to  '  and  the  change  due  to  seepage  to  '  as 
shown  in  Fig.  2.7e  or 

cr  *  =  a  '  +  Acr  '  +  to  '  (2.37) 

f  o  ex  s 

The  stresses  induced  due  to  seepage  may  be  taken 
into  account  by  giving  heed  to  the  water  pressure  distri¬ 
bution.  Consider  the  triangular  element  shown  in  Fig.  2.8. 

Let 

ur  =  resultant  vector  of  forces 


I  . 

- 


' 
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EFFECTIVE  EFFECTIVE 
STRESSES  STRESSES 


FIG.  2.7  SUPERPOSITION  OF  STRESS  STATES  FOR  SOLUTION 
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FIG.  2.8  BOUNDARY  PRESSURES  ON  ELEMENT 
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Ygat  =  saturated  unit  weight 
Ysub  =  submerged  unit  weight 
A  =  area  of  element 
T  =  unit  vector  in  y  direction 

The  change  in  body  force,  AB.F.,  due  to  seepage, 
from  the  initial  condition  of  no  seepage  will  be. 


AB.F .  —  Final  body  force  -  Initial  body  force 

~  ^  A  'feat  ^  +  ur^F  ~  ^-A^sat  +  ur^I 
=  -AYsat  5*  +  “r  +  AYsub  " 

=  "ATW  T  +  ur  (2.38) 


Above  the  phreatic  surface  u  =0,  /.AB.F.  =  -Af  J 

3T  W 

For  no  flow  ur  -  +AYWJ/  /.AB.F.  =  0 

The  AB.F.  found  can  be  equally  divided  among  the 
3  nodes  of  the  element.  AB.F.  is  identical  to  the  seepage 
force  in  the  medium.  The  equations  of  the  change  in  body 
force  are  developed  and  verified  in  comparison  to  the 
seepage  force  in  Appendix  A. 

Both  the  seepage  program  described  earlier  and  the 
stress  program  used  were  left  basically  unchanged.  They 
were  coupled  through  the  use  cf  common  statements  and  several 
additional  subroutines.  One  finite  element  grid  was  used 


for  the  stress  analysis  and  another  was  used  for  the  seepage 
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In  order  to  facilitate  the  assignment  of  body 
forces  from  the  seepage  to  the  stress  grid  and  the  appoint¬ 
ment  Oi  permeabilities  from  the  stress  to  the  seepage  grid/ 
the  element  configurations  were  made  to  correspond  below 
the  phreatic  surface.  Above  this  boundary  the  stress  grid 
also  possessed  the  additional  domain  II  as  shown  in  Fig. 

2 • ?d  with  a  unit  weight  of  yw.  By  having  the  element 
configurations  correspond  exactly  for  both  the  seepage  and 
stress  grids  below  the  free  surface  location  the  delineation 
of  the  phreatic  boundary  was  achieved  in  the  stress  grid. 
Just  as  the  nodes  below  the  free  surface  in  the  seepage 
grid  were  "strung"  along,  both  the  nodes  above  and  below 
this  line  in  the  stress  grid  were  made  to  be  strung  along 
between  the  horizontal  limits  of  the  top  and  bottom  of  the 
excavation.  Fig.  2.9  shows  the  general  method  and  Fig.  2.10 
shows  the  internal  phreatic  surface  location  in  the  stress 
grid  for  one  case  of  the  model  studied.  Since  all  the 
slope  face  nodes  changed  position  due  to  the  moving  free 
surface,  the  release  forces  had  to  be  recalculated  for  each 
stress  analysis.  The  release  force  equations  are  given  in 
Appendix  A. 

To  find  the  correct  free  surface  a  process  consisting 
of  a  double  iteration  technique  was  used.  The  iterations  in 
one  process  are  called  "iterations"  and  in  the  other  are 
called  "cycles".  One  cycle  began  with  a  stress  analysis 
which  included  the  seepage  body  forces  from  the  previous 
cycle.  The  stresses  were  then  used  to  calculate  the 
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FIG.  2.9  METHOD  FOR  INTERNAL  PHREATIC  BOUNDARY 
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FIG.  2.10  INTERNAL  PHREATIC  BOUNDARY  FOR  ONE  CASE  OF  THE  MODEL  STUDIED 
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permeabilities  of  each  element.  The  free  surface  was  found 
with  these  by  iterating.  Convergence  was  taken  when  the 
movement  between  two  cycles  was  less  than  a  specified 
amount.  Seepage  forces  resulting  from  the  analysis  for  a 
homogeneous  medium  were  used  for  the  initial  cycle.  Fig.  2.11 
gives  a  schematic  representation  of  the  processes  and  Fig. 

2.12  gives  the  general  flow  diagram  for  the  computer  program. 

A  listing  of  the  program  is  given  in  Appendix  B. 


STEP  2  STEP  3  STEP 
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FIG.  2.11  SCHEMATIC  PRESENTATION  OF  SOLUTION 


ITERATION 


CALCULATE  THE  HORIZONTAL  AND 
VERTICAL  STRESS  RELEASE  FORCES 


ASSIGN  SEEPAGE  FORCES  RELATED  TO 
SEEPAGE  GRID  TO  THE  MATCHING  NODES  OF 
THE  STRESS  GRID  MODIFYING  THESE 
FORCES  WERE  NECESSARY  FOR  THE 
DISPLACEMENT  BOUNDARY  CONDITIONS 


!  FORM  FINITE  ELEMENT  STRESS  EQUATIONS  i 


SOLVE  EQUATIONS  TO  FIND  AJX,  A Jy 


USE  Cx,  tfy,  A^y  TO  FIND  THE 

AVERAGE  PRINCIPAL  STRESS  FOR 
_ EACH  ELEMENT _ 


FIG.  2.12 


FLOW  DIAGRAM 


CYCLE 


CHAPTER  III 


RELATION  BETWEEN  PERMEABILITY  AND  STRESS 

3 . 1  General 

Until  recently  the  variation  of  permeability  with 
stress  in  rocks  has  been  investigated  mostly  in  the  context 
of  oil  reservoir  engineering.  The  effect  of  all-round 
confining  pressure  has  been  studied  by  Fatt  and  Davies  (1952) 
Knutson  and  Bohor  (1962)  and  others.  The  effect  of  permea¬ 
bility  under  triaxial  conditions  has  also  been  studied  by 
Mordecai  and  Morris  (1971) .  In  other  areas  Brace  et  al 
(1968),  Patching  (1965)  and  Bernaix  (1969)  have  also 
demonstrated  that  permeability  changes  with  stress.  All 
workers  have  found  reductions  in  permeability  as  loading 
increased.  The  fall  was  rapid  at  first  and  changed  less 
rapidly  for  higher  stresses.  In  most  cases  the  reservoir 
sandstone  permeability  reduction  was  less  than  50%.  The 
coals  tested  by  Patching,  the  granites  investigated  by 
Brace  and  the  gneisses  researched  by  Bernaix  showed  much 
larger  reductions.  From  these  studies,  it  was  concluded 
that  this  reduction  occurs  due  to  the  change  in  the  secondary 
structure  of  the  material. 

The  aforementioned  laboratory  determinations  can  be 
supplemented  now  by  permeability  values  from  the  field.  Sher 
man  and  Banks  (1970)  have  summarized  many  known  field  permea¬ 
bility  values  obtained  from  packer  tests.  The  variations  in 
fluid  conductivities  also  arise  due  to  a  secondary  structure 


40 


41 


in  the  rock.  However  this  structure  is  on  a  different 
scale  than  that  in  the  laboratory  specimens. 

In  this  work  it  is  assumed  that  flow  of  water 
through  jointed  and  fissured  rock  obeys  equation  (2.9)  on 
a  macroscopic  scale.  One  intuitively  feels  that  this  is 
correct.  There  may  be  an  analogy  here  with  the  discretization 
of  the  finite  element  approach  where  an  overall  result  is 
correct  but  with  discrepancies  occurring  on  the  element 
scale.  Snow  (1967)  gives  the  justification  for  using 
equations  based  on  Darcy's  law  for  fractured  media.  This 
means  that  the  rock  is  taken  to  conduct  water  everywhere 
and  not  just  at  the  joints  and  fissures.  It  will  now  be 
assumed  that  this  conductivity  at  every  point  is  in  some 
way  related  to  the  effective  principal  stresses  at  every 
point  and  independent  of  shear.  However  it  must  be  remembered 
that  the  permeability  of  interest  arises  due  to  the  secondary 
structure  in  the  rock  and  that  the  primary  permeability  is 
assumed  to  be  small  in  comparison.  This  chapter  attempts 
to  arrive  at  a  reasonable  form  of  equation  relating  permea¬ 
bility  and  stress. 

3 . 2  Development  of  Relation  Between  Permeability  and  Stress 

An  insight  into  the  permeability  of  jointed  rock  may 
be  obtained  using  a  procedure  suggested  by  Terzaghi  (1960) . 

A  model  is  assumed  containing  three  similar  sets  of  joints 
intersecting  each  other  at  right  angles.  The  joint  spacing 
is  D  and  the  width  of  joint  is  b.  From  the  laws  of  laminar 


. 
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flow  of  water  between  parallel  plates  the  average  velocity 
V ,  at  which  water  flows  between  the  plates  is  given  by 


V 


liLbh 

12  y 


(3.1) 


where  y  =  dynamic  viscosity  of  water 


If  f  1  ov/  takes  place  through  the  model  through  two 
sets  of  joints  and  at  right  angles  to  the  third  set  the 
quantity  of  seepage  through  a  square  BxB  is 


Q  =  y?.?.2.  =  i'Ywb^  x  2B^_  _  kj_B2 
D  1 2  y  D 


(3.2) 


where  k 


Ywk^ 

6  y  D 


(3.3) 


Substituting  the  appropriate  values  of  y  and  Yvr ,  k  becomes 


k (cm/sec)  =  16300  b3/D 


(3.4) 


The  relation  between  b,  D  and  k  is  shown  in  Fig.  3.1. 
From  equation  (3.4)  and  the  plot  of  Fig.  3.1,  it  is  observed 
that  if  D  is  constant  the  permeability  in  an  idealized 
model  is  dependent  upon  the  crack  opening  to  the  third 
power.  This  opening  in  turn  will  be  influenced  by  the 
state  of  stress  in  the  region  of  the  crack.  Hence  it  is 
obvious  that  if  the  crack  opening  is  in  any  way  sensitive 
to  stress  the  medium  will  tend  to  have  large  variations 


EQUIVALENT  PERMEABILITY  (cm/sec) 
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FIG.  3.1  RELATION  BETWEEN  JOINT  PARAMETERS  AND  PERMEABILITY 
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in  permeability. 

Fissures  in  a  rock  material  can  be  visualized  as 
discontinuities  with  two  of  the  dimensions  much  greater 
than  the  third.  Walsh  and  Brace  (1966)  have  estimated 
this  aspect  ratio  to  be  in  the  order  of  103  for  low  pressures. 
The  orientation  of  these  planar  openings  may  vary  in  a 
random  manner.  If  most  of  the  head  loss  occurs  due  to  flow 
through  the  fissures  it  is  possible  that  the  permeability 
may  also  be  a  function  of  an  average  crack  width  to  a  power 
although  as  in  the  joint  model  this  too  would  be  an  ideali¬ 
zation  and  oversimplification  to  some  degree. 

In  order  to  arrive  at  some  workable  first  order 
approximation  of  the  permeability  dependence  on  stress, 
again  an  isotropic  jointed  rock  or  an  isotropic  fissured 
rock  material  is  assumed  to  exist  for  the  medium.  As  a 
direct  result  of  this  isotropy  and  the  assumption  of 
conduction  everywhere  a  relation  may  be  established  between 
an  average  principal  effective  stress  and  an  isotropic 
permeability : 


k  —  kx  —  ky  —  f  ( )  (3.5) 

It  may  be  instructive  to  relate  the  average  crack 
width  b,  to  the  average  principal  stress.  The  following  is 
not  presented  as  a  rigourous  formulation  but  only  in  the 
hope  that  it  will  lead  to  a  clearei:  understanding  of 
permeability  dependance  on  stress.  In  general  the  joint  or 


. 
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crack  opening  will  decrease  with  an  increase  in  stress. 
Since  no  data  was  found  for  a  relation  between  stress  and 
crack  opening  a  general  function  could  be  used  to  describe 
the  typical  behaviour.  Fig.  3.2  shows  some  of  the  types 
of  behaviour  that  might  possibly  occur.  It  is  a  type  of 
decay  behaviour  that  a  simple  exponential  function  could 
represent  in  some  cases.  However  the  property  of  an 
applied  exponential  function  is  that  the  rate  at  which  the 
crack  width  changes  is  proportional  to  the  width  itself 
and  this  may  not  be  general  enough.  It  was  found  that  a 
function  of  the  form 


A' 


b 


(3.6) 


where  A'  ,  T  and  M  are  constants 


will  approximate  a  great  variety  of  decay  behaviour.  For 
example,  as  shown  in  Fig.  3.2,  T  =  10  psi  and  M  =  3.0 
could  represent  a  crack  closing  very  quickly  with  stress  or 
T  =  100  psi  and  M  =  1.0  might  represent  a  slowly  closing 
crack  having  dislodged  grain  material  in  the  crack. 

If  equation  (3.6)  is  valid,  equation  (3.3)  can  be 
extended  to  the  form 


k 


A 


(3.7) 


average  (j  +T)N 


This  form  will  hold  true  where  the  main  permeability  arises 
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FIG.  3.2  CRACK  BEHAVIOUR  ACCORDING  TO 
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from  the  secondary  structure  of  the  material  and  alterations 
in  permeability  result  from  the  mechanical  changes  in  crack 
width  due  to  the  applied  stress.  A  relation  that  follows 
this  equation  will  produce  a  straight  line  on  a  log-log 
plot  when  k  is  plotted  against  J±+T.  It  should  be  noted 
that  the  pressure  distribution  is  dependent  on  the  relative 
values  of  permeability  and  hence  will  be  independent  of  A. 

The  quantity  of  flow  is  dependant  on  the  absolute  permea¬ 
bility  and  may  be  scaled  up  and  down  by  changing  A  accordingly. 

3 . 3  Stress  Dependent  Permeability  Data 

The  stress  dependence  of  permeability  in  some  rocks 
has  been  demonstrated  by  various  workers.  The  most  relevant 
data  will  be  reviewed  here.  Eernaix  (1969)  tested  micro- 
fissured  materials  in  an  experimental  set-up  as  shown  in 
Fig.  3.3a  with  typical  results  as  shown  in  Fig.  3,3b. 

Patching  (1965)  tested  the  permeability  of  coal  under  an 
all-round  confining  pressure  with  results  as  shown  in 
Fig.  3.4.  Sherman  and  Banks  (1970)  have  tabulated  data 
which  illustrate  the  decrease  in  permeability  as  measured 
by  packer  tests,  with  depth.  The  data  for  granites  and 
gneisses  are  shown  in  Fig.  3.5. 

In  the  walls  of  the  sample  of  Fig.  3.3a  the  prin¬ 
cipal  effective  stresses  arising  from  fluid  body  forces 
vary  in  some  manner  across  the  sample  wall.  Bernaix  has 
calculated  the  principal  stresses  in  the  sample  assuming 
a  homogeneous  material.  However  in  actuality  these  stresses 
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FIG.  3.3a  TESTING  CONFIGURATION 


FIG.  3.3  VARIATION  IN  PERMEABILITY  WITH  STRESS  (AFTER 
BERNAIX ,  1969) 
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FIG. 


CONFINING  PRESSURE, p.s.i. 


Permeability  vs  confining 
pressure,  samples  from 
one  seam. 


CONFINING  PRESSURE  ,  psi 

FIG  .  3.4b  Permeability  vs  confining  pressure, 

samples  from  several  seams. 


FIG 


3.4 


PERMEABILITY  OF  COAL  (AFTER  PATCHING,  1965) 


DEPTH,  FT 
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.  3.5  DECREASE  IN  COEFFICIENT  OF  PERMEABILITY  OF 
GRANITES  AND  GNEISSES  WITH  DEPTH 
(AFTER  SHERMAN  AND  BANKS ,  1970) 


FIG 
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are  affected  by  the  stress  dependent  permeability  and 
this  makes  it  difficult  to  obtain  a  unique  relation  between 
permeability  and  average  principal  stress.  Strictly  speaking 
the  curves  should  be  viewed  as  an  index  in  evaluating  the 
material.  The  author  has  defined  a  sensitivity  index  as: 


S  -  M~1  bar) 

k  (4-50  bars ) 


(3.8) 


In  Patching's  work  dry  nitrogen  gas  was  used  as 
one  of  the  permeants  in  the  friable  coal.  He  used  "semi¬ 
anthracite  to  low  volatile  bituminous  coaid'.  Field  samples 
were  cut  into  uniform  prisms,  surrounded  by  a  flexible 
epoxy  resin  and  stressed  isotropically  in  an  oil  filled 
pressurized  cell.  No  mention  of  sample  orientation  is 
made.  Permeabilities  were  calculated  from 


yLQ 

(P1-p2)A 


(3.9) 


where  k  is  the  permeability  in  Darcies 

]J  is  the  dynamic  viscosity  in  poises 
L  is  the  length  of  the  specimen  in  cm. 

0  is  the  average  flow  rate  in  cc/sec. 
pl-p2  ps  the  driving  pressure  in  atmospheres 
A  is  the  cross-sectional  area  of  the  sample 
in  cm^ 
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Although  a  gas  was  used,  the  use  of  a  liquid  would 
change  the  viscosity  and  specific  gravity.  Also  gas  slippage 
would  not  occur.  This  latter  phenomena  known  as  the 
"Klinkenburg  effect"  is  a  slipping  of  the  fluid  along 
the  pore  walls  which  gives  rise  to  an  apparent  dependence 
of  permeability  on  pressure.  Mordecai  and  Morris  (1971) 
have  found  this  effect  to  be  of  limited  importance.  In 
any  case  this  results  in  only  a  small  change  in  permeability 
compared  to  changes  caused  by  stress  variation.  Hence 
water  would  give  the  same  type  of  variation  of  permeability 
as  the  gas  used. 

Although  the  permeability  is  somewhat  stress  path 
and  time  dependent  the  overall  behaviour  is  clear.  In 
general  a  curve  can  be  approximated  by  equation  (3.7) 
when  the  initial  segment  at  low  stress  is  curved  with  the 
curvature  decreasing  for  increasing  stress  and  the  curve 
becoming  essentially  linear  for  high  stresses.  It  may  be 
seen  from  Fig.  3,4  that  this  behaviour  is  characteristic 
of  some  of  the  curves  although  all  the  curves  do  not 
follow  this  description.  Plotting  curve  no.  99  a  straight 
line  approximation  may  be  made  for  the  curve  using  T  =  200 
psi  and  giving  N  =  3.36  as  shown  in  Fig.  3.6. 

Sherman  and  Banks  give  a  summary  of  the  known  field 
data  of  the  permeability  of  jointed  rock  with  most  of  it 
being  obtained  from  Snow  (1968a) .  The  permeabilities  were 
determined  by  water  pressure  tests  with  packers.  The 
permeabilities  are  calculated  on  the  assumption  that  the 
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FIG.  3.6 


PATCHING'S  CURVE  NO.  99  FITTED  TO  EQUATION  (3.7) 


' 


flow  is  laminar,  and  that  the  material  is  homogeneous  and 
isotropic  with  respect  to  the  permeabilities.  The  compu¬ 
tation  is  made  according  to: 
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Q  =  k  H  S 

where  k  is  the  permeability  in  cm/sec. 

.  o 

Q  is  the  flow  in  cm  /see. 

H  is  the  net  head  at  the  test  section  in  cm 
S  is  the  shape  factor  in  cm 


The  field 
be  interpreted  by 


data  tabulated  by  Sherman  and  Banks  may 
making  the  following  assumptions: 


a.  The  water  table  is  at  the  surface  (the  net  head, 

H,  is  not  given)  . 

b.  The  horizontal  stresses  =  the  vertical  stresses, 
i.e.  Kq  =  1.0. 

c.  The  tests  are  performed  in  a  flat  lying  area,  i.e., 
the  vertical  stresses  can  be  computed  from  the 
overburden . 

d.  The  packer  tests  have  no  effect  on  permeability. 

e.  The  joint  system  is  isotropic  and  homogeneous. 

f.  The  unit  weight  of  the  rock  yt  =  165  pcf . 

It  is  probable  that  all  of  these  assumptions  will 
be  violated  to  some  degree.  However  they  may  be  used  in 
order  to  evaluate  the  existence  of  a  first  order  approxima¬ 
tion  between  the  average  principal  stress  and  the  permeability. 
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The  results  are  shown  for  items  which  best  fit 
equation  (3.7)  in  Figures  3.7  through  3.20.  From  a  total 
of  22  items  of  3  points  or  more,  greater  than  one  half  fit 
the  equation  to  some  extent,  whereas  the  remaining  data 
was  too  scattered  to  be  used.  In  some  cases  the  curves 
we i- e  fitted  to  a  smooth  fitting  line  passing  through  the 
data  points  for  T  =  0  psi.  In  the  other  cases  the  trial 
value  of  T  was  added  directly  to  the  J1  value  of  the  point. 
It  may  be  observed  that  some  scatter  does  occur  but  for  the 
information  shov/n,  the  data  tends  to  follow  equation  (3.7) 
to  some  degree.  The  relations,  normalized  at  10  psi  so 
uhat  the  permeabilities  in  the  range  of  10  psi  to  an  upper 
limit  may  be  compared,  are  plotted  on  Fig.  3.21.  Fig.  3.22 
shows  T  plotted  against  log  N.  The  occurrence  of  large  T 
usually  results  in  large  N.  Large  values  of  T  will  occur 
when  for  T  =  Opsi  the  curvature  in  the  field  relations 
exist  at  relatively  high  stresses  on  a  log-log  permeability- 
stress  plot. 

3 . 4  Sensitivity  of  Permeability  to  Stress 

In  an  initial  analysis  it  was  noted  that  for  a 
particular  KQ  and  depth  from  the  top  of  slope  to  the  lower 
boundary  a  certain  range  of  effective  stresses  occurred. 

Fig.  3.23  shows  log  permeability  versus  average  principal 
stress  plotted  for  equation  (3.7)  for  several  values  of 
T  and  N  with  the  values  of  permeability  again  normalized 
at  10  psi.  It  may  be  seen  that  the  curve  for  a  certain 
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FIG.  3.7  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS ,  1970) 


PERMEABILITY  (cm/sec) 
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FIG.  3.8  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS ,  1970) 
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FIG.  3.9  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS,  1970) 


PERMEABILITY  (cm/sec) 
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FIG.  3.10 


FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 
(FROM  SHERMAN  AND  BANKS,  1970) 


PERMEABILITY  (cm/sec) 
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AVERAGE  PRINCIPAL  EFFECTIVE  STRESS  (psi) 


FIG.  3.11  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS ,  1970) 
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FIG.  3.12  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS,  1970) 
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AVERAGE  PRINCIPAL  EFFECTIVE  STRESS  (psi) 


FIG.  3.13  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS ,  1970) 
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AVERAGE  PRINCIPAL  EFFECTIVE  STRESS  (psi) 


FIG.  3.14  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS ,  1970) 
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AVERAGE  PRINCIPAL  EFFECTIVE  STRESS  (psi) 


FIG.  3.15  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS ,  1970) 
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FIG.  3.16  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS ,  1970) 


PERMEABILITY  (cm/sec) 
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FIG.  3.17  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS ,  1970) 


PERMEABILITY  (cm/sec) 
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AVERAGE  PRINCIPAL  EFFECTIVE  STRESS  (psi) 


FIG.  3.18  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS,  1970) 
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AVERAGE  PRINCIPAL  EFFECTIVE  STRESS  (psi) 


FIG.  3.19  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS,  1970) 
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FIG.  3.20  FIELD  PERMEABILITY  DATA  FITTED  TO  EQUATION  (3.7) 

(FROM  SHERMAN  AND  BANKS,  1970) 
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FIG.  3.21  COMPARISON  OF  PERMEABILITY  RELATIONS 
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FIG.  3.22  SUMMARY  PLOT  OF  T  AND  N  VALUES 

(FROM  FIELD  DATA  TABULATED  BY  SHERMAN  AND  BANKS, 
1370) 


RELATIVE  PERMEABILITIES  NORMALIZED  AT  10  psi 
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FIG.  3.23  CORRESPONDENCE  FOR  DIFFERENT  PAIRS  OF  T  AND  N 

VALUES 
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combination  of  T  and  N  will  correspond  approximately  to 
another  curve  with  a  different  combination  of  T  and  N  and 
thau  this  correspondence  occurs  in  a  certain  pressure 
range.  For  example,  with  N  =  1.5  and  T  =  20  psi  a  corres¬ 
ponding  variation  occurs  in  the  range  of  10—150  psi  with 
N  —  2.0  and  T  =  40  psi.  It  was  found  while  doing  a  pre¬ 
liminary  investigation  that  the  pressure  distribution 
depended  mostly  on  the  sensitivity  of  the  permeability  to 
stress.  As  a  result  for  the  seepage  analysis  it  was 
decided  to  use  a  combination  of  T  and  N  that  would  give  a 
typical  medium  sensitive  variation  and  a  typical  strong 
variation  of  permeability.  The  curves  chosen  are  shown 
superimposed  as  dotted  lines  on  the  field  data  of  Fig.  3.21. 

3 . 5  Limitations 

With  the  data  at  hand  it  would  appear  that  a  good 
first  approximation  of  permeability  variation  can  be  obtained 
using  equation  (3.7).  It  should  be  kept  in  mind,  however, 
that  a  multitude  of  assumptions  were  used  in  arriving  at 
this  conclusion  and  that  alternative  forms  may  be  equally 
good.  However  this  form  does  fit  much  data  and  is  convenient 
to  use. 

Snow  (1968a)  has  computed  the  minimum  spacing  of 
open  fractures  and  found  that  from  water  pressure  tests 
for  35  damsites  from  the  surface  to  200  ft.  depth,  the 
minimum  spacing  of  open  fracture  increased  from  4  to  14 
ft.  Setting  D  =  Bd+C  in  equation  (3.4),  the  appropriate 
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equation  to  describe  permeability  may  have  the  form 


(J  +T)N  (Bd+C) 
where  d  is  depth 
and  B  and  C  are  constants 

Since  it  was  assumed  that  J  =  f(d,KQ),  it  can  be  seen 
from  the  form  of  equation  (3.11)  that  this  effect  could 
possibly  be  incorporated  in  the  values  of  T  and  N  deter¬ 
mined  previously  and  may  account  for  the  group  of  very 
sensitive  relations  in  Fig.  3.21.  If  a  stress  increase 
occurs  equation  (3.11)  may  be  valid  since  no  new  fractures 
are  formed  and  the  existing  ones  become  tighter.  If  a 
stress  release  occurs  previously  imperceptible  and  closed 
fractures  will  possibly  open  up  so  that  equation  (3.7) 
could  be  valid. 

Although  equation  (3.7)  was  used  for  the  analysis 
the  form  may  be  too  compact.  Besides  using  the  average 
stresses  shear  effects  have  not  been  incorporated  into  the 
equation  mostly  due  to  the  lack  of  data.  For  an  area  of 
high  shear  stresses  where  the  shear  strength  is  approached 
there  would  be  a  tendency  of  movement  along  the  joints  and 
fissures.  A  dilation  effect  would  occur  in  the  area  of 
the  fracture  due  to  its  roughness  increasing  the  average 
width  of  the  crack  resulting  in  increased  permeability. 
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Fig.  3.7  shows  this  effect  in  a  sandstone  (Mordecai  and 
Morris,  1971).  Therefore  the  results  in  Chapter  IV  would 
apply  to  materials  which  have  no  volume  change  with  shear. 


FIG.  3.24  EFFECT  OF  SHEAR  ON  PERMEABILITY 

(AFTER  MORDECAI  AND  MORRIS,  1971) 
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CHAPTER  IV 


RESULTS 

4 . 1  Description  of  Cases  Analyzed 

As  pointed  out  in  Chapter  II  it  was  felt  that  the 
influence  of  stress  dependent  permeability  would  be  most 
pronounced  in  a  steep-walled  excavation.  Hence  only  one 
configuration  of  a  90  slope  was  studied.  The  influence  of 
varying  the  slope  angle  may  also  reveal  important  effects. 

For  a  configuration  considered  to  have  infinite 
base  of  depth  H  and  infinite  crest  of  slope  of  height  H, 
Dunlop ,  Duncan  and  Seed  (1968)  have  shown  that  negligible 
influence  is  felt  in  the  area  of  the  toe  for  an  upslope 
boundary  at  a  distance  of  6H  from  the  toe  and  a  downslope 
boundary  3H  from  the  toe.  For  theoretical  seepage  analyses 
in  slopes,  4  times  the  net  head  or  slope  height  from  the 
toe  to  the  full  upstream  head  is  a  common  boundary  condition 
(Khan,  1971;  Sharp,  1970).  Also  Kealy  and  Busch  (1971) 
using  their  finite  element  seepage  model  have  shown  for  a 
homogeneous  isotropic  embankment  on  a  base  with  a  depth 
approximately  equal  to  the  embankment  height  that  a  base 
to  embankment  permeability  ratio  of  0.1  yielded  almost 
identical  results  to  the  case  when  the  ratio  was  0.001. 
Knowing  that  in  the  present  case,  the  permeability  decreases 
with  depth  in  addition  to  having  the  preceding  guidelines 
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in  mind,  one-half  of  a  centre-line  symmetric  excavation 
as  shown  in  Fig.  4.1  was  chosen  for  the  stress  and  seepage 
analyses.  The  100  ft.  slope  height  was  chosen  arbitrarily. 
This  can  also  be  a  parameter  to  vary  in  future  studies 
since  unfortunately  the  results  of  this  analysis  cannot 
be  normalized  although  they  may  be  extended.  Normalization 
is  prevented  because  the  permeability  variation  occurs  in  a 
particular  stress  range.  However  a  substantial  amount  of 
the  variation  will  occur  for  this  slope  height.  As  a 
result  the  effects  may  be  magnified  in  order  that  the 
governing  principles  in  these  situations  will  be  illustrated 
more  clearly. 

The  stress  grid  shown  has  483  nodes  and  871  elements 
while  the  corresponding  seepage  grid  below  the  free  surface 
has  422  nodes  and  755  elements.  Triangular  elements  were 
utilized  because  a  better  approximation  to  the  permeability 
variation  can  be  obtained  for  the  same  number  of  nodes. 

For  example,  four  nodes  produce  only  one  rectangular  element 
whereas  the  same  number  produces  two  triangular  elements. 

It  is  not  the  number  of  elements  but  the  number  of  nodes 
that  determines  in  large  part  the  computation  time. 

The  grid  was  designed  on  the  basis  of  several 
considerations.  In  addition  to  a  desire  for  detail  at  the 
toe,  an  adequate  representation  of  the  effect  of  the  per¬ 
meability  values  throughout  the  medium  as  discussed  in 
section  2.5  is  required.  The  grid  design  is  also  influenced 
by  a  desire  to  obtain  a  minimum  value  for  a  matrix  quantity 


. 


STRESS  GRIO 


73 


FIG.  4.1  GRID  AND  DIMENSIONS  USED  FOR  ANALYSES 
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known  as  the  "bandwidth"  which  evolves  as  a  characteristic 
feature  of  the  system  of  equations  (2.21)  and  which  affects 
the  storage  requirements  and  computation  time  in  the  finite 
element  analysis. 

It  was  decided  to  vary  the  initial  insitu  stress 
conditions  with  the  variation  in  permeability.  A  simple  K0 
initial  effective  stress  condition  was  chosen  as  shown  in 
Fig.  2.7b.  In  general  however  it  must  be  remembered  that 
finite  horizontal  stresses  often  occur  at  the  surface  for 
rock  materials  (Hast,  1967)  and  this  may  be  another  aspect 
to  investigate  in  the  future.  It  is  common  knowledge  that 
changing  Poisson's  ratio,  v,  in  a  finite  element  analysis 
causes  only  minor  variations  in  the  stress  distribution. 

The  modulus  of  elasticity,  E,  effects  the  displacements 
and  has  little  effect  on  the  stresses  in  plane  strain. 

The  four  cases  studied  and  presented  in  section  4.3 

are : 

1.  A  medium  sensitive  field  permeability  relation, 
i.e.,  T  =  40  psi,  N  =  2.0  and  KQ  =  0.5. 

2.  T  =  40  psi,  N  =  2.0  and  KQ  =  2.0. 

3.  A  sensitive  field  permeability  relation,  i.e., 

T  =  20  psi,  N  =  3.0  and  KQ  =  0.5. 

4.  T  =  20  psi,  N  =  3.0  and  KQ  =  2.0. 

Since  the  pressure  distribution  is  independent  of 
the  factor  A  in  equation  (3.7)  a  value  was  selected  to 
yield  equivalent  permeabilities  for  both  functions  at  0 
psi  average  principal  effective  stress  for  comparison 
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purposes.  A  value  of  10  arbitrary  permeability  units 
(l/t)  at  0  psi  was  chosen.  This  resulted  in  A  =  0.16 
for  T  =  40  psi  and  N  =  2.0  and  A  =  0.8  for  T  =  20  psi 
and  N  —  3.0.  Hence  the  two  permeability  functions  shown  in 
Fig.  3.21  and  used  for  the  analysis  were 

a.  K  =  0.16/(J1+40  )2 

b.  K  =  0.8/  (J,+20  ) 3 

The  seepage  analysis  was  considered  as  a  plane  flow 
problem  and  the  stress  analysis  was  determined  for  plane 
strain  conditions.  The  material  properties  were  taken  as 
modulus  of  elasticity,  E  =  108  psf 
Poisson{s  ratio,  v  =  0.2 

unit  weight  of  rock  material  yT  =  165  psf. 

4.2  Computational  Procedures 

As  outlined  in  section  2.6,  to  initiate  the  solution, 
one  iteration  of  a  seepage  analysis  was  performed  prior  to 
the  first  complete  cycle.  An  initial  position  of  the  free 
surface  was  chosen  which  would  hopefully  tend  to  accelerate 
the  convergence  of  the  subsequent  cycle  process.  Homogeneous 
and  isotropic  permeabilities  were  used  for  the  first  itera¬ 
tion.  The  position  of  the  free  surface  used  for  this 
iteration  was  always  higher  than  the  equilibrium  value 
for  the  homogeneous  and  isotropic  case  and  as  a  result  the 
solution  yielded  negative  pressures  at  and  near  the  free 
surface  nodes.  The  pressures  throughout  the  medium  were 
then  used  to  compute  the  body  forces  for  the  stress  analysis 
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of  the  initial  cycle.  The  results  of  the  stress  analysis 
were  used  in  turn  to  calculate  the  permeabilities.  Hence 
these  initial  permeabilities  were  computed  on  the  basis  of 
an  improper  pressure  distribution  and  a  probable  irregular 
location  of  the  free  surface.  Even  though  this  solution 
initiation  procedure  delayed  cycle  convergence  it  was  a 
necessary  part  of  the  convergence  process.  The  following  cy¬ 
cles  utilized  pressure  distributions  which  had  been  calculated 
in  association  with  free  surfaces  where  little  excess 
negative  or  positive  pressures  usually  existed. 

A  tolerance  value  of  0.1%  of  the  slope  height  or 
net  head  was  generally  used  to  determine  the  free  surfaces 
for  each  cycle.  Since  the  final  position  of  the  free 
surface  was  not  known  the  primary  estimate  was  often  far 
from  an  equilibrium  value.  Hence,  at  times  for  the  initial 
cycles  the  tolerance  value  was  not  met  since  an  inadequate 
number  of  iterations  was  occasionally  specified.  This 
tended  to  result  in  a  requirement  for  additional  cycles 
for  cycle  convergence.  Because  of  initial  difficulties 
due  to  poor  element  shapes  and  because  of  a  lack  of  familia¬ 
rization  with  cycle  convergence  behaviour,  the  specified 
number  of  iterations  was  not  a  constant  and  was  sometimes 
inadequate  as  mentioned  but  in  retrospect  a  value  of  10 
might  have  been  ideal.  However  for  cycle  convergence  the 
tolerances  prescribed  were  always  met  for  at  least  the  last 
two  cycles. 

Usually  the  final  free  surfaces  converged  to  an 
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'-quilibr lum  position  requiring  4  —  8  cyclss.  For  the  initial 
cycles  the  stable  free  surface  for  each  cycle  tended  to 
oscillate  about  the  final  equilibrium  value.  This  oscillating 
behaviour  was  weak  for  most  cases.  However,  strong  oscil¬ 
lating  behaviour  did  occur  for  the  sensitive  permeability 
variation  coupled  with  Kq  =  2.0. 

For  cycle  convergence  a  movement  criteria  was 
used.  If  the  change  between  the  last  two  cycles  performed 
for  any  one  of  the  nodes  on  the  free  surface  was  less  than 
tely  0.1%  of  the  slope  height,  the  final  free 
surface  location  was  accepted.  The  case  described  with 
strong  oscillating  effects  was  accepted  at  0.24%  of  the 
slope  height  after  9  cycles. 

As  a  result  of  the  slow  convergence  in  some  cases, 
an  allowance  was  incorporated  into  the  program  for  performing 
sever al  cycles,  examining  tne  state  of  convergence  and 
continuing  the  converging  process  after  adjusting  the  free 
surface  to  the  value  found  by  the  previous  cycle.  This 
procedure  was  accomplished  in  part  by  reading  and  writing 
the  permeability  values  on  a  tape  file. 

An  IBM  360/67  digital  computer  was  employed  for  the 
calculations.  The  approximate  duration  of  central  processing 
unit  time  for  the  stress  analysis  of  one  cycle  was  50 
seconds  and  for  one  seepage  iteration  was  10  seconds. 

For  each  analysis  the  relevant  data  was  output  on 
punched  cards.  In  addition  to  the  input  data,  element 
permeabilities,  element  stresses  and  nodal  potential  values 
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were  among  the  more  important  data  output.  This  information 
was  used  to  plot  equipotentials ,  permeability  contours  and 
average  principal  effective  stress  contours. 

An  interpolation  and  plotting  program  was  written 
to  reduce  the  output  seepage  data  and  plot  the  equipotentials 
and  the  final  grid  shape.  A  listing  is  given  in  Appendix  C. 

Each  equipotential  plot  results  from  characteristic 
permeability  variations.  Contours  of  equal  permeability 
were  initially  attempted  to  be  output  using  a  library 
contouring  program  but  it  gave  extremely  poor  results  due 
to  an  uneven  density  of  points  within  the  grid  boundary.  As 
a  result  the  equipotential  plot  program  was  modified  to 
compute  the  average  nodal  permeability,  calculate  and  plot 
selected  contours  of  equal  isotropic  permeability.  A 
similar  procedure  was  used  to  output  the  contours  of  the 
average  principal  effective  stresses. 

4 . 3  Results  of  Analyses 

Results  are  presented  for  the  cases  listed  in 
section  4.1  in  Figures  4.3  to  4.6.  For  each  case  the 
equipotential  distribution  represented  by  contours  in 
percentages  of  the  total  head,  selected  contours  of  isotropic 
permeability  and  the  related  contours  of  average  principal 
effective  stress  are  presented.  Fig.  4.2  gives  the  equi¬ 
potentials  in  a  homogeneous  medium  for  comparison.  A 
comparison  of  the  free  surface  locations  are  presented  in 
Fig.  4.7. 
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FIG.  4.2  EQUIPOTENTIALS  IN  A  HOMOGENEOUS  MEDIUM 
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FIG.  4.3a  EQUIPOTENTIALS,  T  =  40  PSI 
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FIG.  4.3b  CONTOURS  OF  ISOTROPIC  PERMEABILITY.  T  =  40  PSI 
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FIG.  4.3c  CONTOURS  OF  AVERAGE  PRINCIPAL  EFFECTIVE  STRESS,  T=40  PSI 
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FIG.  4.4a  EQUIPOTENTIALS ,  T=40  PSI 
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FIG.  4.4c  CONTOURS  OF  AVERAGE  PRINCIPAL  EFFECTIVE  STRESS ,  T=40  PSI 
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FIG.  4.5a  EQUIPOTENTIALS ,  T=20  PSI 
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FIG.  4.5b  CONTOURS  OF  ISOTROPIC  PERMEABILITY ,  T=20  PSI 
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FIG.  4.5c  CONTOURS  OF  AVERAGE  PRINCIPAL  EFFECTIVE  STRESS,  T=20  PSI 
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FIG.  4.6a  EQUIPOTENTIALS ,  T=20  PSI 
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FIG.  4.6b  CONTOURS  OF  ISOTROPIC  PERMEABILITY.  T=20  PSI 
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FIG.  4.6c  CONTOURS  OF  AVERAGE  PRINCIPAL  EFFECTIVE  STRESS,  T=20  PSI 
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FIG.  4.7  COMPARISON  OF  FREE  SURFACE  LOCATIONS 
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To  deduce  the  contour  plots  of  stress  and  permea¬ 
bility,  average  noda.1  values  were  employed.  These  average 
values  may  give  poor  results  under  certain  conditions. 

Each  nodal  value  is  computed  directly  by  taking  an  average 
value  of  the  surrounding  elements.  Cases  where  these 
averages  do  not  give  completely  representative  values  at 
the  nodes  will  occur  when  several  of  the  elements  attached 
to  one  node  physically  extend  more  in  the  direction  of 
maximum  change  of  quantity  contoured  than  the  other  elements 
producing  a  weighted  average.  In  other  situations,  an  irre¬ 
gular  number  of  elements  may  be  attached  to  one  node  due 
to  a  boundary  or  because  of  grid  design  and  may  cause  non- 
representative  values.  At  a  boundary  the  averaging  procedure 
will  be  most  affected  when  the  contour  is  parallel  to  it 
and  will  be  least  affected  when  the  contour  is  orthogonal 
to  it.  For  example,  when  a  contour  approaches  a  boundary 
at  a  low  angle  poor  results  will  occur  yielding  a  type  of 
refraction  effect  close  to  the  boundaries.  In  contrast  to 
these  element  averages,  nodal  potential  values  are  derived 
directly  from  the  finite  element  solution  and  have  not  been 
obtained  through  an  averaging  procedure.  This  will  result 
in  well  defined  values  for  the  equipotentials  and  the  free 
surface . 

As  a  result  of  the  averaging  procedure  somewhat 
improper  or  discontinuous  values  are  obtained  in  the 
following  areas  of  the  contour  plots: 

a.  Within  an  element  dimension  from  both  the  free 
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surface  and  the  upper  grid  boundary  in  the  stress 
plots. 

b.  At  the  lower  boundary  of  the  stress  plots  and  in 
some  cases  in  the  permeability  plots  where  contours 
exist  within  an  element  dimension  from  the  lower 
boundary. 

c.  In  the  area  of  the  toe  and  at  an  approximate 
location  of  a  slope's  height  beneath  the  toe  where 
a  zig-zag  feature  of  the  contour  lines  occurs. 

In  some  of  these  areas  the  proper  value  has  been  drawn  in 
as  a  dotted  line. 

Characteristic  features  may  be  observed  from  the 
equipotential  distributions.  It  is  seen  that  the  phreatic 
surface  is  higher  for  the  more  sensitive  permeability 
function  for  both  KQ  =  0.5  and  KQ  =  2.0.  In  addition  the 
phreatic  surface  is  higher  for  K  =  2.0  than  for  K  ^  =  0.5 
in  both  cases.  The  fact  that  the  permeability  variations 
used  have  more  influence  on  the  free  surface  than  the  KQ 
values  used  may  also  be  noted.  The  equipotentials  are 
more  " S"  shaped  near  the  upstream  boundary  of  the  slope 
for  both  the  more  sensitive  permeability  function  and  the 
higher  KQ  value. 

In  the  contour  plots  of  permeability  and  stress 
distinctive  effects  may  also  be  noted.  Since  permeability 
contours  are  directly  derived  from  the  stress  distribution, 
all  observations  made  for  the  former  will  be  applicable 
to  the  latter.  For  the  higher  value  of  KQ  the  permeability 
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contours  in  the  slope  are  flatter  and  appear  almost  horizontal 
while  for  the  lower  value  the  contours  are  at  more  of  an 
angle  but  they  too  tend  to  become  closer  to  horizontal 
towards  the  back  of  the  slope.  The  permeability  contours 
become  straighter  and  more  horizontal  with  depth  in  ail 
cases  although  poor  average  nodal  values  may  obscure  this 
behaviour.  Into  the  excavation,  away  from  the  face  and 
down  from  the  floor,  the  permeability  contours  become  more 
horizontal  for  all  cases. 

Interesting  results  are  observed  in  the  region  of 
the  toe.  In  all  cases  there  is  an  area  of  low  permeability 
in  this  vicinity.  For  KQ  =  2.0  this  area  extends  a  distance 
directly  beneath  the  toe  while  for  K0  =  0.5  the  area  of 
lower  permeability  extends  down  as  well  but  more  into  the 
slope  than  for  KQ  =  2.0.  This  area  of  low  permeability  is 
less  pronounced  for  KQ  =  2.0  than  for  KQ  =  0.5.  Also  at 
the  toe,  away  from  the  face,  and  below  the  floor  an  average 
principal  stress  decrease  occurs  from  the  original  value 
with  a  greater  effect  occurring  for  the  lower  value  of  KQ . 

With  the  higher  KQ  value,  essentially  a  log  type  of 
decrease  in  permeability  occurs  with  depth.  As  well  as 
this  decrease,  for  KQ  =  0.5  permeabilities  also  become 
lower  when  progressing  in  position  from  inside  the  slope 
to  the  face.  Furthermore,  the  higher  KQ  case  has  an  area 
of  greater  permeabilities  directly  above  the  region  of  low 
permeabilities  at  the  toe  since  a  decrease  in  stress  occurs 
here.  However,  for  KQ  =  0.5  little  of  this  effect  appears 
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to  occur  above  the  toe.  It  is  of  interest  to  note  that  an 
increase  in  curvature  of  the  free  surface  is  observed  in 
this  area  for  KQ  =  2.0  compared  to  KQ  =  0.5.  In  addition 
the  exit  points  approximately  coincide  for  the  sensitive 
relation  and  are  in  close  proximity  to  each  other  for  the 
medium  sensitive  relation. 

.  .  -5 

A  permeability  of  10  x  10  arbitrary  units  should 
occur  at  the  surface  at  the  back  of  the  slope  if  the  average 
principal  stress  is  zero.  The  program  results  reveal  an 
area  of  tension  here  (which  may  be  due  to  the  boundary  con¬ 
dition  at  the  back  of  the  slope)  and  consequently  will 
cause  greater  permeabilities  than  10  x  10" 5  arbitrary  units. 
This  effect  is  more  pronounced  for  the  low  KQ  cases  but  is 
somewhat  obscured  by  poor  nodal  values  at  this  upper  boundary. 

It  is  of  merit  to  observe  that  only  minor  charac¬ 
teristic  differences  may  be  distinguished  in  the  permeability 
contours  between  the  results  of  the  two  sensitivities  used 
for  similar  values  of  KQ. 

In  the  stress  plots,  the  contours  above  the  free 
surface  appear  more  closely  spaced  than  below  it.  This  is 
a  reflection  of  the  total  unit  weight  of  material  acting 
above  the  free  surface  and  the  buoyant  unit  weight  acting 
beneath  it. 
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CHAPTER  V 


DISCUSSION  AND  CONCLUSIONS 

5 . 1  Discussion 

It  is  somewhat  unfortunate  that  the  results  presented 
for  the  iterative  procedure  have  not  or  could  not  be  com¬ 
pared  to  any  existing  solution,  field  situation  or  experi¬ 
mental  setup.  The  solutions  will  have  to  be  accepted  on 
the  merit  of  the  following.  The  basic  approach  of  the 
iterative  method  is  simple  and  the  consequent  results  seem 
entirely  reasonable.  Also  the  finite  element  method  in 
general  has  proven  itself  many  times  and  in  addition  the 
particular  seepage  program  used  was  checked  against  a  close 
form  solution  as  shown  in  Chapter  II. 

The  results  of  section  4.3  illustrate  the  principles 
or  mechanisms  involved  in  the  location  of  the  free  surfaces. 
In  conjunction  with  the  effect  of  the  sensitivity  of  the 
permeability  function,  there  appear  to  be  three  major 
mechanisms  involved  in  raising  the  free  surface  with 
different  combinations  being  accentuated  under  different- 
initial  stress  conditions.  Firstly  there  is  a  continuous 
decrease  in  permeability  with  depth  away  from  the  face. 

In  all  cases  as  the  distance  into  the  slope  from  the  face 
increases  the  permeability  contours  become  increasingly 
more  level  and  in  the  limit  will  be  completely  horizontal 
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with  an  associated  uniform  decrease  in  permeability  with 
depth  due  to  increased  overburden  stresses.  Secondly, 
an  area  at  the  toe  exists  which  has  relatively  lower 
permeabilities  than  the  surrounding  region  due  to  the 
stress  concentration.  This  area  has  a  downward  component 
from  the  toe  under  all  conditions  and  a  horizontal  component 
into  the  slope  under  conditions  of  low  initial  stress. 
Thirdly  a  decrease  in  permeability  occurs  on  a  horizontal 
path  in  the  slope  towards  the  face  as  the  contours  are 
sloped  up  due  to  the  excavation  towards  the  toe  and  face. 

In  cases  of  low  KQ  this  effect  is  greater  than  for  high 
Kq.  All  three  mechanisms  are  dominant  to  some  extent 
with  all  effects  appearing  influential  for  low  K  .  As 
the  initial  stresses  increase  the  first  mentioned  effect 
becomes  more  dominant  than  the  following  two  in  influencing 
the  pressure  distribution.  It  would  appear  that  any  other 
mechanisms  if  they  do  operate  are  of  limited  importance  in 
keeping  up  the  free  surface.  Fig.  5.1  illustrates  the 
mechanisms . 

Stress  release  influences  the  pressure  distribution 
to  varying  degrees  in  two  distinct  areas.  Stress  release 
occurs  away  from  the  toe  and  in  the  excavation  but  this 
has  little  effect  since  most  of  the  head  loss  has  occurred 
prior  to  this  region.  It  does  appear  however  for  cases 
with  high  Kq 1 s  that  stress  release  may  have  some  influence 
on  the  free  surface  location  at  and  towards  the  face.  A 
comparison  of  these  areas  in  Fig.  4.7  shows  for  the  two 
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II.  PERMEABILITY  BLOCK  AT  THE  TOE 


III.  DECREASE  IN  PERMEABILITY  TOWARDS  THE  FACE 


FIG.  5.1  MECHANISMS  WHICH  INFLUENCE  THE  FREE  SURFACE  LOCATION 
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initial  stress  conditions  for  the  same  permeability  function 
that  the  exit  points  are  almost  identical.  The  effect  for 
the  lower  sensitivity  is  less  because  the  entire  free 
surface  and  flow  regime  is  lower  and  farther  removed  from 
this  area  of  stress  release.  The  release  has  resulted  in 
increased  curvature  of  the  free  surface  towards  the  face 
with  resultingly  steeper  gradients  toward  this  vertical 
wall . 

The  major  area  of  permeabilities  in  the  slope  which 
then  may  be  important  in  determining  the  pressure  distribu™ 
tion  appears  to  be  behind  the  face  and  above  a  position 
somewhat  down  from  the  toe.  If  this  is  true  the  variation 
in  permeability  in  this  area  may  be  a  controlling  factor 
and  an  important  variable  in  studying  the  effect  of  stress 
dependent  permeability  in  slopes. 

A  fixed  location  for  the  upstream  boundary  condition 
was  chosen  for  the  study.  However,  climatic  factors  in  the 
hydrological  cycle  as  discussed  in  Chapter  I  will  pro- 
foundly  affect  this  condition.  Terzaghi  (1962)  shows  how 
precipitation  and  runoff  may  affect  the  water  table,  as 
shown  in  Fig.  5.2,  in  close  proximity  to  the  face  due  to 
tension  cracks  at  the  surface.  An  area  of  tension  in  this 
section  of  the  slope  is  most  clearly  illustrated  in  Fig. 

4.6c  in  an  indirect  fashion.  This  additional  mechanism 
may  then  compound  the  effects  of  an  already  high  water 
table  in  undermining  the  stability  of  the  slope. 
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(a)  Location  of  lowest  and  highest 
water-table  in  jointed  rock  with 
uniform  but  low  coefficient  of 
permeability; 

(b)  water-table  in  the  same  rock 
during  heavy  rainstorms  or  snow 
melt  if  permeability  in  space  a  b  c 
is  much  higher  than  elsewhere. 

Line  e  d  represents  the  water-table 
at  times  when  the  joints  drain 
freely,  and  e'd1,  when  they  are 
plugged  with  ice. 


(b) 
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FIG.  5.2 


CLIMATIC-INDUCED  CHANGES  IN  THE  CONFIGURATION 
OF  THE  WATER  TABLE  (AFTER  TERZAGHI ,  1962) 


In  this  work  an  isotropic  linear  elastic  analysis 
has  been  used  and  this  linear  isotropic  property  has  been 
attributed  to  the  rock.  However,  it  is  well  known  that 
non-linear  properties  occur  in  rocks  due  to  the  cracks  and 
joints  (e.g.  Morgenstern  and  Phukan ,  1966;  Zienkiewicz 
and  Stagg ,  1967).  Furthermore,  from  the  crack  behaviour 
shown  in  Fig.  3.2  it  is  evident  that  a  non-linear  stress- 
strain  behaviour  is  also  inherent  in  the  rock  material 
considered . 

It  may  be  argued  that  the  non-linearity  will  not 
nullify  the  general  principles  exhibited  here  since  this 
curvilinear  behaviour  commonly  occurs  in  the  lower  stress 
ranges.  The  area  surrounding  the  toe  being  the  general 
area  of  concern  may  be  considered  to  be  out  of  the  lower 
stress  ranges  with  low  range  effects  in  other  parts  of  the 
slope  having  minimal  influence  in  the  area  of  the  toe. 
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Similarly  it  can  be  argued  for  consistency  that  the  curves 
of  Fig.  3.2  could  be  approximated  by  a  straight  line  at 
higher  stresses.  Since  many  idealizations  have  been  made 
m  this  work  the  inclusion  of  non-linear  behaviour  may  be 
too  much  of  a  sophistication. 

Related  to  the  non-linearity  discussed  but  separate 
from  it  is  a  possible  increase  in  elastic  modulus  with 
depth  (e.g.  for  soft  rocks  Weigh  and  Greenland,  1965). 

Again  as  a  first  approximation  the  modulus  has  been  taken 
to  be  constant. 

Isotropic  rock  properties  have  also  been  assumed 
in  this  work  and  it  may  be  that  anisotropic  effects  are 
important.  However,  Duncan  and  Goodman  (1968)  have  analyzed 
a  jointed  rock  representing  it  by  an  equivalent  orthotropic 
continuum.  They  found  that  the  assumed  orthotropic  nature 
did  not  give  much  different  results  from  where  the  rock 
was  assumed  isotropic.  Hence  the  assumption  of  isotropy 
appears  to  be  very  reasonable. 

The  results  of  this  study  are  valid  essentially  for 
a  slope  height  of  100  ft.  These  results  may  be  extended 
however  by  realizing  that  a  slope  of  a  different  height 
having  a  permeability  relation  that  would  produce  the 
same  permeabilities  as  for  the  100  ft.  slope  would  then  be 
valid.  Thus  the  results  will  hold  for  a  whole  series  of 
permeability  functions  with  the  appropriate  slope  heights. 
For  any  slope  height,  h,  the  results  will  be  valid  for  the 
following  permeability  relations: 
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(Jx+0. 4h) 2 


(5.1a) 
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A 

(Jx+0. 2h) 3 


(5.1b) 


It  is  seen  that  to  arrive  at  the  same  pressure  distribution 
in  higher  slopes  a  less  sensitive  relation  is  required  and 
for  lower  slopes  a  more  sensitive  relation  is  required.  It 
appears  that  the  effect  would  be  less  for  higher  slopes 
with  the  sensitivities  used  and  possibly  more  for  lower 
slopes.  Hence  for  any  sensitivity  there  is  a  slope  height 
where  the  effect  will  be  a  maximum  since  for  lower  slopes 
less  relative  variation  occurs  and  it  is  the  relative 
variations  that  cause  the  effects.  A  plot  of  exit  point 
height  against  slope  height  would  have  the  characteristic 
behaviour  as  shown  in  Fig.  5.3  and  a  family  of  curves 
could  be  obtained  by  varying  either  T  or  N  or  both  in  some 


manner . 


FIG.  5.3  PROPOSED  CHARACTERISTIC  RELATIONSHIP  BETWEEN 
FREE  SURFACE  EXIT  HEIGHT  AND  SLOPE  HEIGHT 
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For  a  higher  slope  with  a  water  table  occurring 
below  the  surface  with  the  same  net  head  and  upstream 
boundary  condition  as  in  this  study  and  with  the  same 
permeability  functions  the  elevating  effect  of  the  free 
surface  would  be  less  since  all  the  stresses  affecting  the 
permeability  will  be  shifted  up  causing  the  variation  to  be 
less . 

The  principles  discussed  initially  will  certainly 
be  in  operation  for  slopes  with  different  geometric  features. 
To  assess  the  effect  of  decreased  slope  angle  from  this 
study  is  difficult,  although  this  may  possibly  be  facilitated 
by  comparing  the  average  principal  effective  stresses  for 
such  cases  with  the  cases  studied.  It  is  likely  that  the 
mechanisms  revealed  here  may  predominate  differently  for 
different  slope  angles.  This  may  also  be  the  case  for 
concave  and  convex  shaped  faces  where  lateral  stresses  will 
tend  to  become  compressive  or  tensile.  Related  to  the 
geometric  configurations  are  the  magnitude  of  the  shearing 
stresses.  In  situations  where  shear  affects  the  permeability 
the  mechanism  of  low  permeability  at  the  toe  may  be  neutra¬ 
lized.  If  the  remaining  mechanisms  are  not  affected  the 
least  influence  due  to  shear  would  be  felt  for  the  higher 
KQ  cases  where  this  mechanism  affects  the  pressure  distri¬ 
bution  the  least. 

It  is  of  interest  to  note  that  the  mechanism  at 
the  toe  may  be  analagous  to  that  controlling  the  emission 
of  gas  in  coal  mines  (Patching,  1970).  Gas  is  generated 
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in  the  coal  due  to  the  maturing  process  of  the  material. 

When  coal  is  mined  underground,  the  resulting  hazards  of 
f ire '  explosion  or  outbursts  and  suffocation  are  in  existence. 
Stress  concentrations  in  the  coal  near  a  face  can  be  due  to 
abutment  loadings  or  due  to  the  shape  of  the  opening. 

These  phenomena  may  reduce  permeability  so  that  high  gas 
pressure  gradients  can  exist  near  the  faces  of  these  areas. 

This  damming  effect  in  both  the  slope  and  the  mine  has 
detrimental  effects  on  safety.  A  finite  element  analysis 
of  the  problem  in  this  particular  area  of  mine  safety  may 
prove  fruitful. 

In  view  of  the  preceding  discussions  it  is  felt 
that  if  these  effects  exist  in  rock  slopes  the  stability 
of  these  slopes  will  be  influenced  to  some  extent  by  this 
behaviour.  However  the  effect  may  often  be  obscured  in 
the  field  due  to  the  natural  non-homogeneity  of  the  material. 

In  some  cases  the  resultant  pressure  distribution  may 
decrease  the  stability  substantially  and  possibly  should 
be  considered  among  the  factors  which  influence  the  stabil¬ 
ity  of  a  slope.  To  describe  the  effect  in  a  quantitative 
manner  seems  impractical  in  view  of  the  typical  non-homogeneity 
that  occurs  in  nature.  Nevertheless  if  such  effects  were  to 
be  taken  into  account  in  a  stability  analysis  it  would  be 
required  to  perform  an  analysis  similar  to  the  one  carried 
out  in  this  work. 

It  is  hoped  that  an  extension  of  the  simplified 
theory  of  section  3.2  to  a  more  realistic  approximation 
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using  the  individual  stresses  could  be  developed.  Different 
joint  distribution  and  properties  should  be  taken  into 
account  (Parsons,  1966;  Snow,  1968b)  in  addition  to  shear 
effects  so  that  anisotropy  as  well  as  non-homogeneity  may 
be  considered.  However  lack  of  insitu  permeability  and 
stress  data  would  hamper  this  development  to  some  extent. 

In  the  event  a  more  advanced  theory  is  evolved,  the  program 
with  its  provision  for  determining  the  stresses  at  any 
orientation  and  facility  for  handling  anisotropy  could  then 
be  easily  extended  to  encompass  increasingly  more  realistic 
cases . 

The  inter-element  configurations  above  and  below 
the  free  surface  may  be  seen  in  Fig.  2.10.  This  one  grid 
was  used  for  lo w  free  surface  cases  as  well  as  high  free 
surface  cases.  Low  phreatic  surfaces  will  cause  an  abun¬ 
dance  of  elements  below  the  free  surface  at  the  toe  and 
will  produce  a  fine  grid  in  this  region.  However,  the 
coarser  grid  above  the  free  surface  at  the  face  may  affect 
the  results  in  the  finer  mesh  below  this  boundary  but  the 
effect,  if  any,  is  minimal.  A  comparison  between  the  average 
principal  stress  contours  of  a  case  with  a  high  free  surface 

and  a  low  free  surface  for  a  particular  K  value  show  the 

o 

same  characteristic  behaviour  in  the  shape  of  the  contours 
just  above  the  toe.  If  effects  do  become  significant 
however,  in  the  present  method  this  could  easily  be  recti¬ 
fied  by  using  a  tailored  mesh  although  the  construction 
of  a  grid  for  this  size  of  problem  can  become  time  consuming. 
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An  alternative  approach  would  be  to  develop  a  method 
involving  a  completely  fixed  stress  grid.  Pressures  at 
the  locations  of  the  nodes  in  the  stress  grid  could  be 
found  interpolating  inside  a  seepage  triangle  after  a 
search  for  this  appropriate  element.  A  special  technique 
would  then  also  have  to  be  developed  for  distributing  forces 
for  the  elements  traversed  by  the  free  surface.  The  body 
forces  would  be  assigned  again  according  to  equation  (2.38). 
The  search  routine  would  be  required  once  more  for  locating 
the  centre  of  a  seepage  element  in  the  appropriate  stress 
element  with  subsequent  interpolation  in  the  stress  element 
being  necessary  if  average  nodal  stresses  were  used.  The 
permeability  could  then  be  assigned  and  the  method  continued 
as  described  in  Chapter  II. 

5 . 2  Conclusions 

The  present  stage  of  development  of  the  computer 
program  was  adequate  to  analyse  the  basic  configuration 
of  excavation  studied.  When  required  the  program  may  be 
easily  extended  to  encompass  shear  and  anisotropy  effects. 

For  the  excavation /the  free  surface  location  is 
influenced  by  the  state  of  initial  ground  stress  and  the 
variation  of  permeability.  With  both  increasing  sensiti¬ 
vity  of  the  permeability  function  and  increasing  values  of 
the  initial  horizontal  stresses  the  phreatic  surface  rises 
above  the  location  for  a  homogeneous  isotropic  medium  for 
the  cases  studied.  The  general  increase  in  elevation  of 
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the  free  surface  location  is  a  result  of  three  basic  types 
of  permeability  variations  which  operate  to  varying  degrees. 
Firstly  a  decrease  in  permeability  occurs  with  depth. 
Secondly  an  area  of  low  permeability  exists  at  the  toe  and 
thirdly  a  decrease  in  permeability  occurs  on  a  horizontal 
path  from  inside  the  slope  towards  the  face. 

The  influence  of  excavation  stress  release  causes 
a  decrease  in  the  average  principal  stress  values  from  the 
original  values.  Two  distinctive  areas  influence  the 
pressure  distribution  to  varying  degrees.  These  areas 
occur  at  a  location  in  the  excavation  away  from  the  toe 
and  at  a  position  above  the  toe. 

If  the  free  surface  behaviour  is  not  wholly  negated 
by  the  assumptions  used  at  arriving  at  it,  the  effects  of 
stress  dependent  permeability  may  in  some  cases  have 

significant  influence  in  decreasing  the  stability  of  rock 
slopes . 


PART  TWO 


CHAPTER  VI 

THE  UNIVERSITY  OF  ALBERTA  MULTISTAGE 
DOUBLE  PACKER  PIEZOMETER 

6 . 1  Introduction 

Part  II  now  deals  with  the  measurement  of  water 
pressures  in  the  field  while  Part  I  has  dealt  with  the 
prediction  of  seepage  behaviour  and  the  associated  pressures. 
Measurements  in  the  field  continue  to  play  an  important 
part  in  estimating  and  monitoring  the  stability  of  soil  and 
rock  slopes. 

Many  types  of  piezometer  designs  are  in  existence 
today.  All  piezometer  systems  in  use  currently  may  be 
broadly  categorized  into  open  standpipe  systems  or  closed 
systems.  In  addition,  closed  systems  may  be  of  hydraulic 
or  diaphragm  design. 

Before  a  general  description  is  given  it  should  be 
noted  that  all  piezometer  systems  must  be  sealed  into 
place.  The  pressure  measured  should  be  comparable  to  the 
actual  pressure  at  the  location  if  no  piezometer  installation 
had  been  made.  The  systems  may  be  sealed  in  a  permeable 
layer  or  a  homogeneous  material  which  could  be  part  of  a 
larger  stratum.  Vaughan  (1969)  has  analyzed  these  two  situa¬ 
tions  for  a  single  piezometer  installation  to  determine  the 
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effects  of  grout  permeability  on  piezometer  measurement. 

Open  standpipe  systems  tend  to  be  simple  in  design. 
They  often  consist  of  a  porous  intake  tip  customarily 
surrounded  by  a  sand  filter.  The  intake  section  is  con¬ 
nected  to  a  standpipe  which  rises  to  the  ground  surface. 

A  change  in  pressure  at  the  tip  will  cause  flow  in  or  out 
of  the  system  until  equilibrium  is  reached  with  the  column 
of  water  in  the  standpipe.  The  measured  head  of  water 
above  the  tip  gives  the  pressure  at  this  point.  For  rela¬ 
tively  pervious  soils  open  standpipe  piezometers  tend  to 
function  adequately.  Casagrande  (1949)  developed  a  system 
incorporating  a  large  filter  intake  with  a  small  diameter 
standpipe.  A  bentonite  soil  was  used  to  seal  the  device 
in  place.  With  soft  soils,  some  designs  can  be  pushed  into 
the  ground  with  the  hope  that  it  is  self-sealing .  The 
Geonor  type  (Bjerrum  et  al,  1965)  operates  in  this  manner. 

Hydraulic  piezometers  commonly  consist  of  a  porous 
tip  connected  to  a  manometer  or  Bourdon  gauge  by  one  or  two 
tubes.  The  flushing  of  air  out  of  the  system  is  facilitated 
by  the  use  of  two  tubes.  Because  of  the  nature  of  the 
system  and  its  use  in  partly  saturated  soil,  air  often 
forms  in  it.  Where  it  is  deaired  a  rapid  response  is 
registered  when  a  change  in  pressure  at  the  tip  occurs.  A 
much  slower  response  will  be  noted  for  air  in  the  system. 

The  pressure  at  the  tip  is  obtained  by  the  addition  of  the 
measurement  at  the  surface  to  the  net  head  of  water  between 
the  tip  and  the  gauge.  This  piezometer  type  is  used  mostly 
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in  embankments  and  shallow  foundations.  Examples  of  these 
closed  hydraulic  systems  are  the  United  States  Bureau  of 
Reclamation  plastic  tip  piezometer  for  embankments  and  the 
Bishop  piezometer  (Bishop  et  al.  ,  1963). 

Diaphragm  piezometers  consist  essentially  of  a 
cylindrical  housing  with  a  diaphragm  at  one  end.  A  porous 
stone  separates  the  diaphragm  from  the  activating  water 
stress.  The  pressure  is  measured  basically  in  two  ways. 

In  one  system  the  diaphragm  acts  as  a  valve.  The  diaphragm 
is  shut  by  external  water  pressure  and  opened  by  internal 
fluid  pressure.  To  take  a  reading  the  internal  pressure 
is  continuously  increased.  A  monitored  quantity  of  either 
pressure  or  volume  measured  through  a  return  tube  will 
reveal  when  the  valve  opens.  At  this  point  the  pressures 
on  both  sides  of  the  diaphragm  are  balanced  and  the  activating 
pressure  can  be  read  at  the  surface.  Among  others,  air 
activated  diaphragm  systems  have  been  used  (Warlam,  1963). 

With  the  other  method  an  electrical  conductor  is  used  to 
communicate  the  amount  of  deflection  of  the  diaphragm  to 
the  surface.  The  electrical  currents  can  be  converted  to 
pressure  by  a  previously  determined  calibration  curve.  To 
sense  the  pressure  on  the  diaphragm  electrical  piezometers 
may  have  vibrating  wire  mechanisms  (Scott  and  Kilgour , 

1967)  or  strain  gauge  mechanisms  (Lundgren,  1966). 

Upon  installation  of  electrical  instruments, 
measurements  may  often  be  conducted  in  a  simple  and  efficient 
manner.  Vibrating  wire  devices  operate  on  the  principle 


c 


117 


that  the  natural  frequency  of  a  stretched  wire  is  dependent 
on  the  tension  in  the  wire.  To  obtain  a  reading  a  stretched 
wire  attached  to  the  diaphragm  in  the  instrument  is  activated 
by  an  electro-magnet.  Sensors  in  the  instrument  then  pick 
up  the  frequency.  An  electronic  system  at  the  surface 
monitors  the  frequency  which  is  consequently  related  to 
pressure  by  a  calibration  curve.  Transducer  piezometers 
commonly  use  electrical  strain  indicators  to  gauge  the 
pressure  on  the  diaphragm.  Upon  deflection  of  the  diaphragm 
an  electrical  resistor  mounted  on  it  changes  its  electrical 
conductivity.  'The  resistances  of  the  transducer  and 
strain  indicator  combine  in  circuit  so  that  a  Wheatstone 
bridge  is  developed.  As  a  result  the  strain  indicator 
resistance  reading  can  be  converted  to  pressure  by  a 
calibration  curve. 

One  important  factor  influencing  the  acquisition  of 
pore  pressure  data  is  a  measurement  time  lag  (Hvorslev, 

1951) .  This  time  lag  is  dependent  on  the  permeability  of 
the  soil,  the  shape  and  size  of  the  intake  and  the  volume 
factor  of  the  system,  i.e.,  the  volume  of  water  which  must 
flow  to  or  from  the  system  for  a  unit  change  in  pressure. 
Where  quick  response  times  are  required  for  rapid  determina¬ 
tion  of  pressure  or  where  the  low  permeability  of  a  soil 
may  cause  excessively  long  response  times  it  is  desirable 
to  have  small  volume  factors.  In  general  diaphragm  piezo¬ 
meters  have  this  characteristic.  The  most  convenient 
diaphragm  system  from  the  standpoint  of  measurement  is  an 
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electrically  operated  piezometer.  Since  the  costs  of  these 
electrical  units  in  the  past  have  been  somewhat  excessive 
as  well  as  the  fact  that  diaphragm  systems  tend  to  be 
unreliable  for  very  long  periods  of  time  (Lindberg,  1965), 
many  of  these  types  have  been  developed  to  be  retrievable 
(e.g.  Brooker  et  al,  1968). 

A  recoverable  transducer  device  has  been  designed 
as  shown  in  Fig.  6.1a  which  is  installed  directly  into  a 
drilled  hole  (Lundgren,  1966).  It  is  sealed  by  means  of 
an  inflatible  rubber  bag.  The  bag  is  filled  with  water  by 
admitting  water  through  line  B  as  shown  on  Fig.  6.1a  and 
allowing  air  to  escape  through  line  A.  When  the  water  in 
the  bag  is  placed  under  pressure  the  piezometer  will  be 
sealed  under  proper  conditions  at  the  desired  elevation. 

A  comparison  between  the  response  times  of  this  device  and 
a  porous-tip  piezometer  with  a  null  indicator  is  given  by 
Lundgren  and  shown  in  Fig.  6.1b.  The  two  piezometer 
systems  were  installed  in  a  clayey  soil  at  a  depth  of 
16.5  ft.  The  water  table  was  originally  5  ft.  from  the 
surface.  A  third  borehole  had  been  drilled  so  that  each 
of  the  three  holes  formed  the  corner  of  an  equilateral 
triangle  with  2  ft.  sides.  The  effect  measured  resulted 
from  lowering  the  water  level  10.5  ft.  in  22  minutes.  The 
results  reveal  an  apparently  sensitive  system. 

A  device  similar  to  the  one  described  previously 
was  developed  by  the  University  of  Alberta  in  co-operation 
with  the  Geological  Survey  of  Canada  (Brooker  et  al,  1968). 
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Recoverable  electric  pieiometer  has 
been  developed  for  use  in  the  field,  to  provide 
greater  accuracy,  easier  installation,  and  more 
rapid  response. 


FIG.  6.1a  BOREHOLE  PIEZOMETER 
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Comparative  sensitivity  of  electrical 
transducar  versus  porous-tip  piezometer  and 
null  indicator  is  shown.  Effect  measured  was 
lowering  of  water  level  in  an  adjacent  open 
borehole  10.5  ft  in  22  min. 


FIG.  6.1b  RESPONSE  BEHAVIOR 


FIG. 


6.1  BOREHOLE  PIEZOMETER  AND  RESPONSE  CHARACTERISTICS 
(AFTER  LUNDGREN,  1966) 


12  0 


Known  as  the  UA-GSC  transducer  piezometer,  it  is  shown  in 
Fig.  6.2a-i  and  is  intended  as  an  instrument  for  measuring 
pressures  in  impermeable  clay  shales.  In  this  method  a 
4-6  in.  hole  is  drilled.  A  well  point  is  attached  to  the 
end  of  a  1 h  in.  I.D.  galvanized  iron  pipe  of  borehole 
length  and  lowered  into  the  hole  to  the  required  depth.  A 
sand  filter  is  installed  around  the  wellpoint  and  a  cement 
or  clay  grout  is  placed  immediately  above  the  filter  sand 
to  provide  the  required  seal.  This  installation  in  the 
borehole  is  shown  in  Fig.  6.2a-ii.  The  piezometer  is  lowered 
to  the  well  point  and  the  packers  are  inflated  using  an 
air  reservoir.  The  instrument  will  be  sealed  into  the 
desired  area  when  the  rubber  packers  expand  against  the 
walls  of  the  iron  pipe.  A  comparison  of  response  curves 
between  the  device  and  an  open  standpipe  is  given  by  Brooker 
et  al  and  is  shown  in  Fig.  6.2b. 

Another  system,  ingenious  in  design  and  somewhat 
different  from  the  previous  two  models  has  been  fabricated 
for  use  in  rock  (Londe,  1971).  The  transducer  instrument 
is  known  as  a  continuous  borehole  piezometer.  A  borehole 
is  permanently  lined  with  a  continuous  flexible  rubber 
membrane  as  shown  in  Fig.  6.3a,  closed  at  the  bottom  and 
fitted  at  the  top  with  a  pressure  lock.  The  membrane  is 
filled  with  water  so  that  the  pressure  at  all  points  in  the 
membrane  will  be  greater  than  the  surrounding  pressure  in 
the  rock.  This  results  in  a  seal  against  the  rock  for  the 
entire  length  of  the  membrane.  A  measuring  probe  shown  in 
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FIG.  6.2a  UA-GSC  PIEZOMETER  SYSTEM 


PIEZOMETER  U-65-1 


FIG.  6.2b  RESPONSE  BEHAVIOUR 
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6.2  UA-GSC  PIEZOMETER  AND  RESPONSE  CHARACTERISTICS 
(AFTER  BROOKER  ET  AL,  1968) 
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FIG.  6.3a  PIEZOMETER  SYSTEM 
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FIG.  6.3b  RESPONSE  BEHAVIOUR 


FIG 


6.3  CONTINUOUS  BOREHOLE  PIEZOMETER  AND  RESPONSE 
CHARACTERISTICS  (AFTER  LONDE ,  1971) 
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Fig.  6.3a  carrying  two  packers  to  isolate  a  short  section 
of  the  hole  is  passed  through  the  pressure  lock  and  lowered 
to  any  desired  elevation.  After  activation  of  the  packers 
an  electrically  operated  piston  in  the  probe  can  be  engaged 
so  that  the  pressure  P  inside  the  membrane  at  the  probe 
will  be  less  than  the  pressure  p  in  the  rock.  Water  will 
flow  from  the  surrounding  rock  causing  equalization  of  the 
pressures  on  both  sides  of  the  membrane.  A  pressure  trans¬ 
ducer  in  the  probe  will  measure  the  pressure  in  the  hole. 
Since  the  monitored  pressure  will  tend  towards  the  true 
value  asymtotically ,  the  piston  is  activated  to  bracket  the 
magnitude  of  the  correct  pressure.  Fig.  6.3b  shows  a 
response  curve  using  this  bracketing  technique  in  quartzite. 
With  this  technique  a  continuous  water  pressure  profile  of 
the  borehole  can  be  developed.  The  present  instrument  was 
purported  capable  of  being  used  in  86  and  116  mm  boreholes 
at  depths  of  "several  hundred  meters". 

6 . 2  Description  of  the  University  of  Alberta  Multistage 
Piezometer  System 

A  comparable  device  to  the  continuous  borehole 
piezometer  has  been  developed  at  the  University  of  Alberta 
for  use  in  soil  and  rock  slopes.  The  instrument  can  be 
described  as  a  borehole  double  packer  transducer  piezometer 
and  is  basically  the  UA-GSC  piezometer  with  an  additional 
packer  below  the  transducer.  The  piezometer  unit  is  shown 
in  Fig.  6.4.  The  design  has  evolved  through  the  efforts  of 
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FIG.  6.4  ASSEMBLY  DRAWING  OF  THE  UNIVERSITY  OF  ALBERTA 
DOUBLE  PACKER  TRANSDUCER  PIEZOMETER 


125 


staff  and  students  in  the  geotechnical  section  of  the 
civil  engineering  department  at  the  university.  Piezometer 
components  are  shown  in  Appendix  D. 

The  multistage  system  can  take  readings  at  pre¬ 
selected  levels.  From  the  stratigraphy  revealed  by  the 
borehole  log  of  a  6-8  in.  drill  hole,  the  desired  locations 
of  well  points  in  the  borehole  are  determined.  The  well 
points  are  assembled  with  sections  of  lh  in.  I.D.  polyvinal 
chloride  (PVC)  pipe  so  that  the  measuring  points  will  be 
at  the  proper  elevations  when  the  assembly  is  lowered  in 
the  hole.  The  well  points  consist  of  1  ft.  sections  of 
lh  in.  I.D.  galvanized  iron  casing,  threaded  at  both 
ends,  perforated  throughout  with  approximately  h  in.  holes 
and  surrounded  with  2  layers  of  fine  mesh  wire.  Impermeable 
seals  are  placed  above  and  below  each  wellpoint.  Filter  sand 
is  positioned  at  the  points  and  filler  sand  is  deposited 
between  them. 

Readings  can  be  taken  after  the  installation  is 
completed.  The  piezometer  is  lowered  down  the  PVC  casing 
to  the  desired  wellpoint.  The  flexible  rubber  membranes 
are  inflated  with  compressed  air  to  seal  the  device  at  the 
desired  area  or  layer  developing  in  the  formation  of  a 
closed  measuring  system.  A  diagrammatic  representation  of 
the  piezometer  in  operation  is  given  in  Fig.  6.5.  Readings 
of  pressure  versus  time  will  reveal  the  equilibrium  pressure 
at  the  point. 
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The  effects  of  drilling  mud  are  important  in 
maintaining  the  original  flow  conditions  in  permeable  water 
bearing  strata  after  the  installation  is  completed  and  no 
measurements  are  being  taken.  The  use  of  a  bentonite 
drilling  slurry  during  the  initial  boring  tends  to  penetrate 
and  seal  the  strata  so  that  longitudinal  communication 
between  the  levels  chosen  is  minimized.  The  lew  volume 
factor  of  the  instrument  develops  rapid  response  times  in 
areas  where  the  sealing  mud  has  effected  low  permeability 
in  the  permeable  layers.  Fig.  6.6  shows  an  example  of  a 
possible  equilibrium  flow  regime  arising  in  a  hypothetical 
situation  to  illustrate  the  desired  effects  of  drilling 
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FIG.  6.6 


HYPOTHETICAL  SITUATION  TO  ILLUSTRATE  THE  DESIRED 
EFFECTS  OF  BENTONITIC  DRILLING  MUD 


CHAPTER  VII 


FIELD  INVESTIGATION 

7 . 1  Location  of  Installation 

Investigators  working  in  areas  associated  with  the 
Upper  Cretaceous  Edmonton  Formation  have  found  that  infor¬ 
mation  concerning  possible  or  existing  perched  water  tables 
is  important.  In  a  report  on  a  slope  stability  study  in 
the  immediate  vicinity  of  the  University  of  Alberta  campus, 
E.W.  Brooker  and  Associates  (1970)  express  the  limitations 
of  their  analysis  by  remarking,  "Perched  groundwater  con¬ 
ditions  applied  to  failure  surfaces  extending  deep  into  the 
(Cretaceous)  bedrock  can  lead  to  an  overly  pessimistic 
indication  of  stability  conditions".  In  another  case, 
Eigenbrod  and  Morgenstern  (1971)  analyzing  a  slide  associated 
with  Cretaceous  bedrock  and  perched  water  conditions  show 
that  if  the  water  pressure  distribution  in  the  slide  had 
been  deduced  naively,,  an  improper  interpretation  of  the 
shear  strength  determined  from  the  analysis  would  have  been 
made.  Hence  in  the  past  the  knowledge  of  irregular  water 
conditions  associated  with  Cretaceous  bedrock  has  been 
deemed  important. 

An  appropriate  location  with  Cretaceous  bedrock  was 
sought  to  ascertain  the  working  characteristics  of  the 
piezometer  system,  its  strengths,  its  weaknesses  and  if 
possible  a  reasonable  indication  of  the  groundwater  conditions. 
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Since  no  existing  open  boreholes  were  available  an  attempt 
was  made  to  determine  a  convenient  and  suitable  location. 

In  particular  it  was  decided  to  select  a  site  where  the 
stratigraphy  was  generally  known  and  that  entertained  a 
high  probability  of  having  one  or  more  perched  water  tables. 

As  shown  in  Fig.  7.1  a  location  near  the  University  of  Alberta 
campus  was  chosen  which  would  penetrate  the  Upper  Cretaceous 
bedrock.  A  maximum  depth  of  200  ft.  of  hole  was  planned. 

The  geology  of  the  area  has  been  studied  previously 
and  the  following  brief  description  is  taken  from  the  work 
of  Thomson  (1970).  Deposits  of  four  ages  are  present  in 
this  area.  Starting  with  the  oldest,  these  are,  (1)  Upper 
Cretaceous,  (2)  Saskatchewan  Sands  and  Gravels,  (3)  Glacial 
(dense  till)  and  (4)  Postglacial  or  recent  (lacustrine 
deposits) .  The  strata  of  the  Upper  Cretaceous  consist  of 
interbedded  grey  shale,  grey  sandstone,  bentonitic  sand¬ 
stone,  bentonitic  shale,  coal  and  bentonite.  Approximately 
20  ft.  above  the  river  a  coal  seam  2  to  3  ft .  thick  outcrops. 

A  section  through  the  location  is  given  in  Fig.  7.2.  It  is 

well  known  that  the  coal  beds  tend  to  be  laterally  extensive. 

/ 

Since  the  coal  layers  often  drain  into  the  valley,  irregular 
pressure  distributions  abound. 

7 . 2  Field  Procedures 

A  Failing  1500  drilling  unit  was  used  for  the  work 
associated  with  the  installation.  The  stratigraphy  in  the 
Upper  Cretaceous  bedrock  was  determined  by  obtaining 


4 


. 


2 

O 

H 

Eh 

< 


131 


FIG.  7.1  LOCATION  OF  INSTALLATION  SUPERIMPOSED  ON  STUDY  AREA  OF  THOMSON  (1970) 
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FIG.  7.2  SECTION  THROUGH  LOCATION  FOR  INSTALLATION  (AFTER  THOMSON,  1970) 
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continuous  core  using  both  a  pitcher  sampler  and  a  4  in. 
core  barrel.  Seventy  feet  of  hole  which  was  partly  in  the 
sands  and  gravels  were  cased  with  a  9h  in.  steel  casing 
since  previous  experience  in  these  materials  pointed  to  the 
possible  collapse  of  the  hole  and  to  sampling  difficulties 
caused  by  gravel  falling  to  the  base  of  the  boring.  Drilling 
was  ceased  at  a  depth  of  197  ft.  After  the  lower  127  ft. 
were  reamed  to  8  in. ,  bentonite  drilling  mud  was  circulated 
in  the  hole  until  fluid  losses  ceased.  Clear  water  was 
then  circulated  to  clean  the  boring. 

The  log  of  the  test  hole  is  shown  on  the  right  side 
of  Fig.  7.3.  A  possible  coal  layer  exists  at  a  depth  of 
121  ft.  but  core  recovery  was  not  obtained.  It  was  decided 
to  establish  well  points  for  the  installation  as  shown  on 
the  left  of  Fig.  7.3  at  four  elevations: 

a.  Well  point  no.  1  in  the  Saskatchewan  Sands  and 
Gravels  at  59.  ft.  (Elev.  2177.  ft.) 

b.  Well  point  no.  2  in  the  clay  shale  at  105.  ft. 

(Elev.  2129.  ft.) 

c.  Well  point  no.  3  in  a  4  ft.  coal  stratum  found  at 
156. ft.  (Elev.  2080.  ft.)  This  is  15-20  ft.  above 
the  river. 

d.  Well  point  no.  4  in  a  6-12  in.  coal  layer  found  at 
187.  ft.  (Elev.  2049.  ft.)  This  is  10-15  ft.  below 
the  river. 

The  sealing  material  was  selected  on  the  basis  of 
insitu  soil  conditions,  ease  of  placement  and  permeability 
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characteristics.  A  seal  adequate  for  many  insitu  soil 
conditions  was  desirable.  For  the  depths  associated  with 
the  installation ,  it  was  essential  that  a  pumpable  grout 
be  used  since  tamping  a  seal  into  place  at  these  depths  was 
not  feasible.  Past  experience  and  laboratory  tests  have 
indicated  that  pure  cement  grouts  or  bentonite  cement  grouts 
have  high  permeabilities  relative  to  bentonite  clay.  This 
was  thought  to  arise  in  large  part  from  cracking  due  to 
snrinkage .  As  a  result  it  was  decided  to  use  a  low  permea¬ 
bility  AM-9  grout  (Lamb,  1960)  or  a  cement  grout  with  an 
additive  to  induce  expansion. 

Because  of  the  highly  toxic  nature  of  AM- 9  grout 
and  because  experience  with  its  use  was  lacking,  a  cement 
grout  with  an  expansive  additive  was  used  for  sealing  the 
well  points  in  place.  An  economical  additive  known  as 
'Entraplast  A1*  was  mixed  with  normal  Portland  cement  in  a 
ratio  by  weight  of  1  part  additive  to  100  parts  cement. 

This  grout  has  been  tested  and  used  successfully  by  Lyell 
(1969)  in  a  harder  rock.  For  pure  cement  paste  tested 
under  lab  conditions  where  the  water-cement  (W/C)  ratio  of 
the  paste  was  less  than  0.7  and  0.5,  the  permeability  was 
found  to  be  less  than  10  and  10  cm/sec  respectively 
(Powers  et  al,  1955).  However  increases  in  permeability  in 
the  cement  often  occur  due  to  preferential  flow  paths 
caused  during  settlement  (or  "bleeding'’ )  and  as  mentioned 
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previously  microcracking  associated  with  shrinkage. 

The  pumping  system  on  the  drilling  unit  was  used 
for  placing  the  grout.  The  pump  was  a  duplex  reciprocating 
type  with  possible  pumping  rates  of  up  to  180  gallons  per 
minute.  Since  the  lowest  possible  rate  of  pumping  was  much 
too  high  to  place  the  grout  without  dispersion,  the  pumping 
system  was  modified  with  a  by-pass  valve  to  obtain  minute 
flow  rates  of  grout  into  the  hole.  The  initial  seal  was 
placed  without  a  bypass  system  resulting  most  likely  in  a 
poor  seal.  Approximately  three  seals  were  placed  per  day. 
Fig.  7.4  shows  the  pumping  system  used. 

Grout  was  pumped  down  the  hole  using  a  1  in.  I.D. 

PVC  pipe  in  20  ft.  sections.  Additional  10  ft.,  5  ft., 

2  ft.  and  1  ft.  sections  were  used  to  obtain  any  length  of 
grout  pipe  required.  The  grouting  pipe  was  also  used  as  a 
gauge  for  measuring  the  sand  fill  depth  and  the  grouted 
length  of  seal.  The  depth  measurements  were  facilitated 
by  the  sedimentation  of  a  small  layer  of  pea  gravel  to  form 
a  base  from  which  measurements  could  be  taken. 

Clean  uniform  Ottawa  sand  was  sedimented  to  form 
the  filter  material  at  the  well  points.  A  well  graded  sand 
was  used  as  a  filler  material  between  well  points.  For  the 
present  installation  the  use  of  well  graded  sand  by  sedimen¬ 
tation  caused  segregation  and  a  longer  wait  time  for  settle¬ 
ment  than  was  necessary.  Hence  if  sand  fill  is  required 
in  any  future  installations  a  well  sorted  coarse  material 
should  be  used.  Sand  settling  time  for  200  ft.  of  hole  was 
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such  that  if  installations  deeper  than  this  are  attempted 
a  new  method  of  placing  the  sand  should  be  found  in  order 
to  decrease  stand-by  time  for  machinery  when  high  priced 
rental  equipment  is  used. 

It  was  decided  that  a  set  cement  was  not  required 
for  continuation  of  installation  of  the  seals,  filter  sand 
and  filler  sand.  This  decision  resulted  from  the  considera¬ 
tion  of  several  factors.  Laboratory  tests  revealed  that  a 
minimum  of  "bleeding"  occurred  for  W/C  ratios  of  less  than 
approximately  0.5.  As  a  result,  for  larger  ratios,  the 
effective  stresses  in  the  fresh  cement  will  be  low  since 

relatively  high  pore  water  pressures  in  the  cement  slurry 

/ 

are  induced  due  to  consolidation  and  sedimentation.  Hence 
sand  penetration  for  a  slurry  with  a  W/C  ratio  greater  than 
0.5  will  be  large  due  to  its  inability  to  support  the  sand 
grains  because  of  low  inter-cement  particle  stresses.  It 
was  also  observed  in  the  laboratory  that  the  sand  penetrated 
the  cement  a  small  distance  so  that  a  sand  cement  layer  was 
formed.  This  layer  acted  as  a  filter  squeezing  water  from 
the  cement,  consolidating  it  and  giving  it  strength.  Be¬ 
cause  of  the  confining  action  of  the  sand-cement  layer, 
shear  strength  in  the  grout  was  not  required  to  prevent 
penetration.  The  penetration  pressure  of  the  sand  will 
also  be  reduced  since  the  sand  and  seals  were  to  be  placed 
under  water  so  that  the  buoyant  weights  of  the  sand  and 
cements  are  only  acting.  In  addition  if  sand  is  placed 
over  the  grout,  arching  in  the  sand  against  the  sides  of 


♦ 


B  •'  ^ 


139 


the  borehole  limit  penetration  stresses  to  a  large  extent. 
Hence  with  these  mechanisms  in  operation  sand  penetration 
into  the  fresh  grout  will  not  be  excessive. 

The  following  two  procedures  were  consequently 
observed  in  the  installation  of  the  materials  in  the  field. 
After  grouting,  a  token  wait  time  of  30  minutes  was  allowed 
before  sand  was  sedimented,  so  that  the  cement  had  some 
time  to  settle  if  it  was  dispersed  at  the  top  of  the  seal. 
Also  a  thin  layer  of  fine  sand  was  deposited  between  sand 
and  cement  layers.  This  was  necessary  because  the  coarser 
fraction  in  the  sand  used  settled  first  disturbing  the  seal 
if  a  fine  sand  supporting  base  was  not  present. 

In  addition  to  the  two  preceding  procedures  a  W/C 
ratio  of  0.5  or  less  was  required  for  placing  the  seals. 

The  pumping  system  used  had  a  minimum  circulating  volume  of 
20-25  gallons.  With  this  initial  volume  of  water  a  minimum 
of  six  80  lb.  bags  of  cement  had  to  be  mixed  at  any  one 
time  to  obtain  a  W/C  ratio  of  0.5.  Cement  was  added  to 
the  circulating  grout  in  the  pumping  system  until  this 
ratio  was  reached.  Water-cement  ratios  were  determined 
from  grout  samples  which  were  taken  and  weighed  in  the  field. 
The  W/C  ratio  was  interpreted  according  to  Fig.  7.5. 

The  grout  volumes  pumped  were  calculated  to  give  a 
5-8  ft.  seal.  Due  to  the  expansive  nature  of  the  grout  it 
was  desirable  to  use  a  seal  as  long  as  feasible  so  that 
the  middle  expanded  against  the  upper  and  lower  parts  of 
the  set  cement.  Seal  lengths  were  gauged  by  monitoring  the 
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level  in  the  mixing  drum.  The  estimated  final  locations  of 
the  seals  are  summarized  on  the  left  of  Fig.  7.3.  For  well 
point  no.  1  the  hole  had  sloughed  in  to  give  a  much  longer 
seal  than  anticipated  and  thus  grouting  in  the  well  screen. 
This  seal  could  have  been  removed  and  replaced  quickly  with 
a  more  flexible  arrangement  of  pumping  system.  However 
since  the  drilling  program  was  behind  schedule  a  lack  of 
time  coupled  with  the  thought  that  the  piezometer  could 
quickly  measure  through  the  grout  resulted  in  no  remedial 
action.  When  dealing  with  perched  conditions  it  is  always 
desirable  to  seal  the  entire  length  of  a  relatively  homo¬ 
geneous  and  impermeable  stratum  if  it  is  of  considerable 
extent.  Hence  if  the  water  pressure  distribution  at  the 
well  point  in  the  shale  is  far  from  hydrostatic  an  incorrect 
reading  may  be  obtained  since  the  hole  was  not  sealed  for 
the  entire  length.  Grouting  was  not  performed  further  on 
both  sides  of  the  point  as  shown  in  Fig.  7.3  because  ini¬ 
tially  a  5-8  ft.  seal  was  planned  for  4  well  points.  As  a 
result  a  lack  of  'Entraplast  A'  was  on  hand  to  grout  the 
entire  length.  This  and  other  deficiencies  in  foresight 
then  resulted  in  the  length  of  seal  at  this  location. 
Nevertheless  interesting  results  from  the  piezometer 
readings  are  obtained  for  well  point  no.  2.  As  shown  in 
Fig.  7.3,  the  location  of  the  bottom  seal  is  higher  than 
it  should  have  been.  This  was  a  result  of  a  mis judgment  in 
the  installation  procedure  of  the  sand.  Since  the  permea¬ 
bility  of  the  coal  layer  is  much  larger  than  the  clay  shale 
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the  position  of  the  bottom  seal  is  not  critical. 

7 . 3  Readings  Obtained 

At  the  outset ,  difficulties  were  encountered  in 
obtaining  reliable  measurements.  The  problems  met,  essen¬ 
tially  fell  into  two  main  categories.  Firstly  a  class  of 
problems  appeared  which  concerned  lowering  the  instrument 
in  the  casing  and  thereafter  finding  the  correct  level  for 
measurement.  A  second  group  of  difficulties  was  related  to 
the  maintenance  of  packer  pressure  in  the  piezometer.  A 
majority  of  the  problems  were  overcome  to  obtain  the 
readings  presented  below. 

Even  before  the  well  points  were  grouted  in  place 
some  difficulties  were  encountered  in  lowering  the  piezometer 
for  the  length  of  the  PVC  casing.  It  was  found  that  the 
instrument  would  not  pass  the  first  coupling.  The  upper 
section  of  pipe  was  dismantled  and  it  was  observed  that  the 
PVC  pipe  ends  had  been  pressed  inward  to  form  a  circular 
protruding  lip.  This  might  have  been  due  to  a  possible 
over-tightening  of  the  coupling.  After  the  edges  were 
filed  flush  the  pipe  was  reassembled  and  the  instrument  was 
lowered  to  the  bottom  without  serious  difficulty. 

Notwithstanding  the  precautions  described  previously, 
during  subsequent  measurements  similar  difficulties  were 
encountered  at  some  locations.  This  was  due  possibly  to 
the  obstruction  described,  to  a  recessed  annular  ring 
formed  at  the  joints  on  which  the  piezometer  could  catch 
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or  to  some  other  block  mechanism.  The  piezometer  always 
passed  the  restricted  locations  but  often  in  the  attempt 
the  sleeve  securing  the  packer  membranes  tended  to  be 
loosened  so  that  the  packers  were  unable  to  maintain 
pressure.  If  pressure  loss  was  observed  the  piezometer 
was  hauled  to  the  surface  for  inspection. 

Usually  when  a  pressure  loss  occurred  the  membrane 
was  immediately  suspected.  On  one  occasion  however  after 
both  membranes  had  been  replaced  and  firmly  secured  a  per¬ 
sistent  pressure  loss  still  occurred.  After  the  instrument 
was  brought  to  the  laboratory  for  inspection  a  close 
examination  revealed  a  leak  in  the  air  line  several  inches 
above  the  connection  of  the  hose  to  the  piezometer.  A 
leak  in  the  pressurizing  hose  was  not  considered  likely 
because  the  hydraulic  line  was  of  double  stainless  steel 
jacketed  design  with  a  nylon  covering.  After  this  section 
of  hose  was  replaced  pressure  loss  problems  were  of  a  minor 
nature. 

The  air  and  the  electrical  lines  weighed  0.152  and 
0.149  lb.  per  ft.  respectively  combining  to  give  0.3  lb. 
per  ft.  In  addition  the  piezometer  body  weighed  approximately 
7  lbs.  As  a  result  for  the  lowest  well  point  approximately 
a  60  lb.  force  was  required  to  hold  the  assembly  in  place 
and  handling  became  awkward  for  this  depth.  Also  it  was 
found  while  marking  off  the  well  point  elevations  on  the 
air  hose  in  the  lab  that  it  stretched  substantially.  Since 
the  electrical  cable  elongated  easily  as  well,  the  ability 
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to  locate  the  well  point  levels  precisely  was  hindered. 

On  one  occasion  a  membrane  was  damaged  when  it  expanded 
into  the  perforations  on  the  inside  of  the  well  point. 

Because  of  this  stretching,  the  piezometer  level  had  to 
be  adjusted  in  small  increments  to  find  the  well  point  level. 
Also  to  aid  in  the  seating  the  distance  between  the  packers 
was  increased  by  incorporating  an  extending  section  in  the 
body  of  the  piezometer.  Since  the  instrument  might  not 
have  been  centred  properly  at  the  well  point  a  maximum  error 
in  measurement  in  the  order  of  1  ft.  of  water  could  have 
occurred  due  to  misalignment.  On  the  average  however  the 
error  was  likely  smaller.  After  the  levels  were  found 
both  the  air  and  electrical  lines  could  be  marked  to  give 
the  proper  elevation  for  future  readings. 

Another  difficulty  arose  when  searching  for  the 
correct  level.  The  compressed  air  reservoirs  with  volumes 
of  288  cu.  in.  could  be  pressurized  to  only  90  psi  with 
existing  facilities  instead  of  an  allowable  400  psi.  Con¬ 
sequently  at  the  lower  well  point  when  packer  pressures  had 
to  be  in  the  order  of  60  psi  or  more  only  one  or  two  infla¬ 
tions  were  possible.  An  attempt  was  made  to  overcome  this 
inconvenience  by  the  use  of  a  highly  pressurized  C02  con¬ 
tainer.  Since  measurements  at  the  time  were  taken  in  sub¬ 
zero  weather  the  temperature  of  the  CC>2  gas  as  it  vaporized 
into  the  system  was  very  low.  In  the  installation  the 
temperature  was  probably  in  the  order  of  40°F  and  subse¬ 
quently  when  the  gas  warmed  a  pressure  increase  occurred. 


- 


145 


This  had  an  undesirable  effect  on  the  measured  response 
curves.  In  addition,  a  clumsy  system  of  valves  on  the 
CC>2  container  resulted  in  a  sudden  high  pressure  on  one 
occasion  rupturing  the  membranes.  With  this  experience, 
the  use  of  CC>2  as  a  pressure  source  was  discontinued. 

A  faulty  regulator  on  one  of  the  air  reservoirs 
was  also  responsible  for  pressure  loss  several  times. 

On  account  of  the  problems  described  previously 
an  initial  set  of  mostly  incomplete  readings  were  acquired. 
These  were  taken  in  December  of  1971.  The  temperature  in 
the  field  at  this  time  ranged  to  -40°F  and  resulted  in  poor 
performance  of  the  equipment  on  the  surface  and  very 
unpleasant  working  conditions.  A  second  set  taken  in  the 
beginning  of  1972  proved  mere  fruitful.  The  level  of  water 
in  the  casing  dropped  between  these  two  times  to  an  equi¬ 
librium  level  12  ft.  above  well  point  no.  3.  When  readings 
were  desired  for  the  upper  two  well  points,  the  casing  had 
to  be  filled  above  the  level  where  the  readings  were  desired. 

Figures  7.7  through  7.10  show  the  measurements 
obtained  at  the  four  well  points  when  interpreted,  according 
to  the  calibration  curve  of  Fig.  7.6.  Generally,  packer 
pressures  15-20  psi  in  excess  of  the  head  above  the  well 
point  were  used  to  seal  the  piezometer  in  place.  As  dis¬ 
cussed,  difficulty  was  experienced  in  obtaining  the  readings 
so  that  all  the  curves  are  not  complete.  Since  great 
efforts  were  made  to  obtain  all  of  these  readings,  they 
will  nevertheless  be  presented. 
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FIG.  7.8  RESPONSE  AT  WELL  POINT  NO 
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RESPONSE  AT  WELL  POINT  NO.  4 


151 


The  curves  obtained  for  well  point  no . 1  are  shown 
in  Fig.  7.7.  The  initial  measurement  after  completion  of 
the  installation  appears  to  show  very  little  head  of  water 
at  the  point.  It  is  interesting  to  note  the  slow  initial 
response  for  the  first  measurement  since  this  part  of  the 
curve  is  almost  horizontal.  The  following  measurement  shows 
a  head  of  34.8  ft.  at  one  period  in  time.  Although  the 
initial  reading  is  not  complete  it  appears  that  the  response 
times  could  be  similar .  Water  had  to  be  added  to  the  casing 
for  the  second  measurement. 

Fig.  7.8  shows  the  curves  procured  for  well  point 
no.  2  in  the  shale.  The  initial  curve  shows  that  the  pos¬ 
sible  equilibrium  value  could  have  been  in  the  order  of 
20  ft.  or  less  with  response  times  in  the  order  of  10,000 
to  20,000  minutes.  The  subsequent  measurement  shows  a 
minimum  at  300  minutes  followed  by  an  increase  in  pressure 
with  increases  still  occurring  at  40,000  minutes  to  a  value 
in  the  order  of  60-70  ft.  of  water.  It  is  important  to 
note  that  before  this  measurement  was  taken  the  well  point 
was  above  the  water  level  in  the  casing.  The  initial  response 
time  appears  much  slower  than  the  following  response. 

The  behavior  for  the  4  ft.  coal  layer  at  well  point 
no.  3  is  shown  in  Fig.  7.9.  Equilibrium  pressures  for  three 
trials  are  generally  reached  at  approximately  2000  minutes 
yielding  a  value  of  11  to  13  ft.  of  head.  An  equilibrium 
value  of  11.0'  was  attained  on  Dec.  13/71  and  a  value  of 
12.4'  on  Feb.  1/72.  It  is  of  interest  to  note  the  much 
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more  rapid  response  for  the  initial  trial  than  for  any  of 
the  other  trials. 

The  equilibrium  pressure  for  well  point  no.  4  shows 
17.7'  of  head  in  Fig.  7.10.  The  6-12  in.  coal  layer  is 
estimated  to  occur  10-15  ft.  below  river  level.  Equili¬ 
brium  occurs  at  approximately  1000  minutes  for  the  second 
measurement.  The  reading  of  Dec.  9/71  was  obtained  soon 
after  the  completion  of  the  installation  and  at  the  time  a 
large  head  of  water  was  still  in  the  casing.  This  is  the 
reason  for  the  large  pressure  dissipation.  For  the  two 
readings  shown  the  response  time  for  the  earlier  reading 
appears  to  be  longer. 
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CHAPTER  VIII 


DISCUSSION  AND  CONCLUSIONS 

8 . 1  Discussion 

With  the  installation  and  operation  of  the  multi¬ 
stage  piezometer  system  insight  has  been  obtained  concerning 
grouting  procedures,  quality  of  the  cement  seals,  piezometer 
operating  characteristics  and  the  pressure  distribution  in 
the  slope.  With  this  knowledge  it  is  felt  that  decisions 
may  be  formed  regarding  future  implementation  of  the  system 
or  a  modified  version  thereof. 

The  grouting  procedure  utilized  may  be  improved  for 
more  rapid  placement  of  seals.  The  cost  of  the  installation 
is  in  part  determined  by  the  placement  rate  when  expensive 
rental  equipment  is  used  and  it  is  desirable  to  keep  all 
costs  minimal.  There  were  two  main  problems  in  positioning 
the  seals.  The  lack  of  knowledge  of  the  seal  height  placed 
and  the  absence  of  a  facility  to  easily  flush  out  seals 
improperly  installed  are  the  two  difficulties  that  ideally 
should  be  overcome  in  future  grouting  operations.  The 
following  discusses  features  which  may  be  incorporated  in 
such  a  system. 

With  the  setup  used,  only  one  pump  was  in  operation. 
When  it  was  necessary  to  circulate  water  through  the  hole 
during  grouting  operations  the  pump  had  to  be  cleared  of 
slurry  and  replaced  with  water.  This  switchover  procedure 
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resulted  in  time  delays  and  caused  dilution  of  the  existing 
grout  when  cement  was  again  circulated.  Hence  it  is  desi¬ 
rable  to  have  an  auxiliary  unit  for  pumping  water. 

A  controlled  or  accurately  determined  grout 
height  is  required  in  the  case  of  a  lower  seal.  A  method 
should  be  included  for  clearing  out  the  grout  pipe  so  that 
the  volume  of  seal  applied  is  known  more  accurately.  For 
example,  in  200  ft.  of  1  in.  I.D.  pipe  approximately  1  cu. 
ft.  of  material  exists.  This  is  a  substantial  amount  and 
is  comparable  to  the  volume  placed  for  the  seal.  If  the 
pipe  is  disconnected  from  the  pump  and  raised  before  it 
is  cleared  out  the  grout  will  add  to  the  seal  and  also 
become  dispersed  with  the  water  in  the  hole.  If  pumped 
water  is  used  in  the  clearing  process  disturbance  and 
dilution  of  the  seal  may  occur.  Hence  compressed  air  might 
be  ideal  for  the  procedure. 

Even  though  attempts  can  be  made  to  accurately 
control  the  seal  length  various  factors  may  cause  error. 
With  the  installation  of  seals  one  is  always  "working  in 
the  dark"  as  it  were  and  hence  it  is  desirable  to  incor¬ 
porate  a  system  to  measure  the  height  of  grouting  applied. 
This  may  possibly  be  performed  using  resistivity  measure¬ 
ments  with  probes  being  attached  to  the  grouting  rod  in 
some  manner.  Thus  having  the  upper  location  under  control 
more  grout  could  quickly  be  added  if  desired  or  the  entire 
seal  could  be  flushed  away  if  the  length  was  longer  than 
predicted  and  might  interfere  with  the  operation  of  the 
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well  point.  Also  dispersion  of  the  sealing  material 
might  be  detected  earlier  if  the  grout  were  pumped  in  place 
too  quickly. 

In  addition  to  installing  the  seals  the  pipe  could 
also  aid  with  the  placement  of  the  filter  and  filler  sands. 
Since  sedimented  sand  is  usually  in  a  loose  state  it  would 
be  advantageous  to  have  it  compacted.  Vibration  of  the  sand 
would  assist  in  this.  Hence  the  grouting  pipe  could  be  used 
to  densify  the  material  by  attaching  a  vibrating  source  to 
it.  It  is  not  desirable  to  vibrate  the  PVC  piezometer 
casing  since  this  may  disturb  the  previous  seals  which  have 
just  begun  to  set. 

In  these  ways  the  grouting  pipe  could  serve  a  multi¬ 
tude  of  functions.  Besides  placing  the  cement  seals  and 
measuring  the  height  of  sand  fill,  it  could  be  used  to 
flush  out  the  hole  quickly  when  required,  compact  the 
sedimented  sand  and  measure  the  length  of  seal  placed.  A 
representation  of  the  proposed  system  is  given  in  Fig.  8.1. 

It  is  very  disappointing  to  note  that  response  times 
measured  with  the  multistage  system  are  much  longer  than 
those  indicated  by  Brooker  et  al  (1968)  in  Fig.  6.2b. 

Brooker  et  al  infer  that  the  single  packer  piezometer  has 
a  very  small  volume  factor  with  resultant  rapid  response 
times.  It  is  true  that  the  volume  factor  of  a  transducer 

_  r  _n 

alone  is  estimated  to  be  in  the  order  of  10  D  to  10 
cm5/'gm  which  is  almost  negligible  compared  to  that  of  a 
closed  hydraulic  system  such  as  a  Bishop  piezometer  system. 
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FIG.  8.1  PROPOSED  MODIFIED  GROUTING  SYSTEM 
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However  the  flexibility  of  the  air  filled  packers  must  also 
be  taken  into  account.  This  will  substantially  increase 
the  amount  of  water  required  to  activate  the  system.  The 
volume  factor  of  the  double  packer  piezometer  was  measured 
using  the  apparatus  as  shown  on  Fig.  8.2.  Care  was  taken 
to  exclude  any  air  that  would  affect  the  results.  From  the 
curves  on  the  figure  it  may  be  seen  that  an  average  volume 
of  0.1  cc  is  required  for  a  pressure  change  of  1  psi  giving 
an  average  volume  factor  of  1.42  x  10“3  cm5/gm  which  is 
comparable  to  the  stiffness  of  a  Bishop  piezometer  system 
with  approximately  200  ft.  of  polythene  leads  (Penman,  1961). 

As  a  result  of  the  information  of  Fig.  8.2  it  can 
be  concluded  with  a  large  degree  of  certainty  that  the 
response  curve  given  by  Brook  et  al  as  shown  in  Fig.  6.2b 
for  well  points  in  shale  is  actually  a  result  of  dissipa¬ 
tion  of  excess  pressure  through  the  cement  seals.  These 
seals  must  have  been  full  of  microcracks  with  resultant 
high  permeabilities.  The  quick  response  is  then  not  a 
result  of  the  purported  small  volume  factor  in  conjunction 
with  the  low  permeability  of  the  shale  but  of  a  pervious 
seal . 

From  the  obtained  response  curves  in  general  and 
in  view  of  the  preceding  it  is  felt  that  the  expansive 
additive  'Entraplast  A'  is  wholly  adequate  for  sealing 
piezometers  into  boreholes.  Its  use  is  recommended  where 
insitu  field  conditions  will  allow  a  rigid  material  to 
provide  adequate  sealing  qualities. 
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VOLUME  CHANGE  (CC)  CORRECTED  FOR  SYSTEM  STIFFNESS  TO  POINT 
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FIG.  8.2  VOLUME  FACTOR  OF  DOUBLE  PACKER  PIEZOMETER 
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From  Fig.  8.2  it  may  be  seen  that  pressure  increases 
will  not  occur  beyond  42  psi  which  indicates  that  for  the 
testing  configuration  used, 8  psi  is  required  to  seal  the 
piezometer  since  the  packer  pressure  was  50  psi.  Hence  this 
was  taken  as  an  indication  of  the  critical  packer  pressure 
to  be  used  in  the  field.  In  the  work  it  was  assumed  that 
on  the  basis  of  the  information  of  Fig.  8.2  which  was 
determined  during  the  field  measurements  that  an  excess 
pressure  in  the  order  of  20-30  psi  would  be  enough  to  seal 
the  piezometer  in  all  cases.  The  curves  obtained  seem  to 
indicate  adequate  seals.  For  example,  curves  for  well 
point  no.  3  were  acquired  under  a  variety  of  initial  heads 
and  appear  to  proceed  generally  to  similar  equilibrium 
values.  Since  criteria  were  not  established  prior  to  the 
field  measurements,  a  fuller  investigation  is  recommended 
to  determine  a  complete  set  of  criteria  in  cases  of  different 
combinations  of  pressures  above,  below  and  between  the 
packers  so  that  field  conditions  may  be  simulated  and  to 
evaluate  other  factors  which  may  possibly  have  an  influence 
in  sealing  the  piezometer  such  as  increasing  the  cable 
tension  after  inflation.  Although  the  author  has  confidence 
in  the  results  obtained,  if  it  is  found  that  the  packer 
pressures  used  were  not  sufficient  in  general,  readings  may 
be  repeated  with  the  appropriate  pressures. 

Response  times  for  piezometer  systems  have  been 
studied  and  analyzed  (Hvorslev,  1951;  Gibson,  1963).  Hvor- 
slev  predicted  for  an  incompressible  material  that  times  to 
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a  certain  percentage  response  will  be  a  constant  for  any 
change  in  pressure  at  the  piezometer  tip.  Possibly  this 
effect  is  brought  out  in  the  response  curves  of  well  point 
no.  3  in  the  coal.  Although  percentage  responses  will  not 
be  the  same  for  the  entire  measurements,  it  is  interesting 
to  note  the  curves  tend  to  come  to  equilibrium  at  a  similar 
time.  Readings  at  the  other  well  points  were  too  scarce 
or  too  varied  to  make  any  such  observation. 

Information  may  be  determined  from  the  equilibrium 
values  of  the  response  curves.  For  well  point  no.  1, 
assuming  an  initial  equilibrium  reading  much  less  than  the 
following  measurement,  it  appears  that  freezing  of  the 
slope  face  may  have  effected  an  area  of  low  permeability 
causing  the  water  table  to  rise  substantially.  The  melting 
effect  of  the  sun  may  have  some  influence  on  determining 
the  upper  water  distribution  in  the  slope.  The  slope  in 
which  the  multistage  piezometer  is  situated  is  on  the  south 
side  and  removed  from  the  effects  of  the  sun.  If  this 
hypothesis  is  valid  the  general  influence  of  freeze-up 
could  possibly  be  determined  by  attempting  a  correlation 
of  slide  frequency  with  the  sun-side  or  shade-side  of  the 
river . 

The  low  permeability  of  the  seal  which  has  grouted 
in  the  screen  of  well  point  no.  1  will  cause  considerable 
time  lag.  Any  change  in  groundwater  conditions  that  do 
occur  will  be  reflected  with  a  large  phase  change. 

Although  the  initial  reading  at  well  point  no.  2 


r 


161 


is  incomplete  as  well,  it  is  interesting  to  speculate  that 
the  freeze-up  also  effects  the  pressure  in  this  area.  It 
is  possible  that  an  installation  time  lag  caused  by  stress 
release  could  have  caused  the  low  initial  reading.  However, 
low  pressures  occurred  simultaneously  at  the  upper  two 
well  points  and  a  water  table  increase  is  then  apparently 
reflected  at  both  points  as  well.  A  large  measurement 
time  lag  appears  at  the  second  well  point  because  the  pres¬ 
sure  became  zero  as  the  equilibrium  water  level  in  the 
casing  fell  below  this  point.  The  well  point  acted  as  a 
drain  with  the  result  that  the  pressure  in  this  area 
became  low.  The  soil  subsequently  must  have  consolidated 
but  to  a  very  small  degree  since  the  shales  are  generally 
quite  hard.  When  the  piezometer  was  installed  the  second 
time,  the  draining  was  stopped  and  water  flowed  into  the 
area  again  increasing  the  pressure.  This  process  may  be 
observed  in  Fig.  7.8  as  it  reached  a  minimum.  Equilibrium 
was  not  yet  attained  even  after  one  month.  The  initial 
response  curve  hopefully  precludes  a  theory  of  a  leaky 
cement  seal.  The  phenomena  described  emphasizes  one  of 
the  weaknesses  of  the  present  multistage  system.  Generally 
an  increased  measurement  time  lag  will  be  effected  in 
impervious  and  compressible  material  when  the  pressure  due 
to  the  column  of  water  in  the  casing  at  the  well  point  is 
lower  or  higher  than  the  pressure  in  the  soil.  A  single 
packer  piezometer  in  impervious  material  may  be  more  prac¬ 
tical  since  a  large  change  in  pressure  at  the  point  relative 
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to  the  pressure  in  the  soil  should  not  occur  as  the  water 
in  the  casing  should  not  drain  out.  Hence  it  is  recommended 
■chat  measuring  points  not  be  placed  in  impervious  materials 
when  used  in  conjunction  with  an  open  multistage  system 

f 

since  measurement  time  lag  will  be  large  even  if  the  volume 
factor  is  small.  As  pointed  out  in  section  7.2,  depending 
on  the  pressure  distribution  in  the  shale,  longer  seals 
might  have  resulted  in  different  equilibrium  readings. 
However,  measurement  time  lag  would  also  have  occurred  for 
that  case. 

The  well  point  in  the  4  ft.  coal  layer  might  also 
bring  out  some  of  the  same  effects  of  freezing  as  may  have 
been  observed  with  the  upper  two  well  points.  According 
to  the  section  of  Thomson  (1970) ,  the  coal  layer  outcrops 
approximately  700  ft.  from  the  installation  so  that  this 
appears  to  be  a  well  draining  layer. 

The  rapid  response  of  the  initial  trial  at  this 
well  point  is  interesting.  The  following  is  a  possible 
explanation.  With  the  initial  high  head  the  drilling  mud 
penetrated  to  a  certain  extent.  As  the  piezometer  was 
installed  and  packers  inflated  the  pressure  decreased  as 
shown  on  Fig.  7.9  and  the  effective  stresses  in  the  coal 
layer  at  the  point  increased.  As  a  result  the  drilling 
mud  was  squeezed  further  into  the  coal  fractures  with  the 
tightening  of  the  cracks.  Apparently  the  effect  of  the 
combination  of  the  drilling  mud  and  the  crack  tightening 
has  affected  subsequent  measurements.  The  fact  that  some 
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fractured  coals  may  have  large  changes  in  permeability  with 
stress  has  already  been  brought  out  in  Fig.  3.4.  One  would 
ordinarily  expect  and  hopefully  obtain  response  characteris— 
tics  as  shown  with  the  initial  measurement. 

The  bottom  coal  layer  appears  to  drain  freely  into 
the  river.  The  layer  is  approximately  10-15  ft.  below  the 
river  level  with  a  measured  head  of  water  of  17.4  ft. 

There  is  likely  flow  out  of  the  4  ft.  layer  and  into  the 
6~12  in.  layer.  The  flow  is  enough  to  account  for  the 
head  above  well  point  no.  3  and  flow  into  well  point  no.  4. 
It  is  possible  that  the  pressure  in  the  lower  layer  may  be 
a  function  of  river  level.  This  could  be  ascertained  within 
the  limits  of  the  instrument  by  continuous  measurements  of 
pressure  and  river  level. 

It  is  recommended  that  measurements  be  continued 
on  the  installation  for  the  acquisition  of  additional  data 
accompanying  seasonal  pressure  fluctuations  as  well  as  to 
verify  the  measurements  already  obtained. 

During  the  drilling  program  a  bentonite  bed  was  not 
observed  above  the  4  ft.  coal  layer  and  may  have  been 
missed.  In  the  work  of  Thomson  (1970)  a  6  in.  bentonite 
layer  is  purportedly  above  this  coal  stratum.  This  may 
act  as  an  impermeable  barrier  depending  on  the  overall 
relative  permeabilities  of  the  bentonite  and  shale.  Eigen- 
brod  and  Morgenstern  (1971)  also  report  a  bentonite  bed 
above  one  of  the  coal  layers.  Although  geological  reasons 
may  exist  for  a  systematical  sequence  of  this  nature. 
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bentonite  seams  will  not  necessarily  overlie  coal  beds  and 
furthermore  could  be  below  the  coal  as  well. 

A  possible  pressure  diagram  for  winter  conditions 
may  exist  as  shown  in  Fig.  8.3.  The  effect  of  a  possible 
coal  seam  at  121  ft.  and  the  consequences  of  possible  inade¬ 
quate  seals  in  the  shale  may  cause  the  distribution  to  be 
considerably  different  but  on  the  basis  of  the  obtained  data 
the  interpretation  is  as  adequate  as  any.  Also  if  bentonite 
acts  as  an  impermeable  barrier,  the  pressure  distribution 
may  be  as  indicated  by  the  dash-dotted  lines  and  could 
appear  as  the  dashed  line  if  the  bentonite  has  no  effect. 

In  terms  of  slope  stability  this  perched  distribution  will 
influence  analyses  to  give  higher  factors  of  safety  if 
failure  surfaces  extend  below  the  perched  water  tables. 

It  is  apparent  that  the  multistage  piezometer 
system  is  not  as  practical  as  was  first  thought.  The  major 
drawback  of  the  present  system  is  the  much  longer  than 
expected  response  times.  This  was  exactly  the  one  aspect 
that  made  the  instrument  appear  practical.  However  the 
system  may  yet  be  useful  for  measuring  pressures  in  perme¬ 
able  strata.  If  compression  of  air  in  the  packers  is  the 
main  cause  of  the  magnitude  of  the  volume  factor  the  incor¬ 
poration  of  an  incompressible  fluid  in  some  manner  may 
decrease  this  value.  For  example,  Lundgren  (1966)  uses 
water  for  inflation  of  the  packers  and  this  may  have 
decreased  the  volume  factor  of  the  instrument  used.  Pre¬ 
ferably  however,  another  packer  system  should  be  designed 
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PRESSURE  IN  FT. OF  WATER 


FIG.  8.3 


POSSIBLE  PRESSURE  DISTRIBUTIONS  IN  SLOPE 
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to  establish  a  low  compliance  of  the  instrument  and  yet 
be  able  to  produce  adequate  seals.  If  such  a  system  is 
developed,  almost  instantaneous  response  will  be  attained 
and  will  make  the  system  immeasurably  more  useful. 

Other  problems  have  also  detracted  from  the  practical 
aspects  of  the  system.  For  depths  greater  than  200  ft. 
problems  will  occur  in  the  lowering  of  the  instrument  due 
to  the  additional  weight  of  the  cables.  If  the  present 
cables  are  used  difficulties  will  also  be  encountered 
finding  the  exact  well  point  location  due  to  the  elasticity 
of  these  lines.  Both  of  these  problems  might  be  overcome 
with  the  use  of  a  relatively  rigid  steel  cable  in  conjunction 
with  a  hoisting  system  to  place  the  instrument.  The  air 
line  could  then  be  replaced  by  a  lighter  rubber  hose. 

In  view  of  the  preceding  other  systems  should  be 
considered.  As  of  late  very  economical  Japanese  built 
pressure  transducers  have  come  on  the  market.  These  could 
be  fitted  into  a  proper  housing  tip  and  permanently  grouted 
into  the  desired  layers.  Although  the  present  devices 
contain  some  hysteresis  and  a  small  amount  of  creep  as 
revealed  by  laboratory  tests,  for  the  pressures  measured 
in  the  study  these  would  have  been  ideal.  However  since 
one  is  always  wary  of  burying  a  wholly  electrical  system 
because  of  their  past  unreliability  over  long  periods  of 
time,  a  long-term  program  should  be  initiated  to  reveal  the 
true  limitations  of  these  transducers  in  the  field,  modi¬ 
fying  them  as  required.  A  system  for  placing  these  could 
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be  designed  as  shown  in  Fig.  8.4.  The  improved  grouting 
procedure  discussed  could  be  used  for  installing  the  seals. 
With  this  piezometer  system  there  will  be  negligible 
measurement  time  lags  in  addition  to  the  convenient  and 
rapid  acquisition  of  data  at  all  points  simultaneously. 

8 . 2  Conclusions 

The  grouting  system  used  was  generally  adequate  to 
install  the  cement  seals.  Some  inconveniences  and  limita¬ 
tions  associated  with  the  technique  have  been  noted  and  a 
modified  method  has  been  proposed. 

The  grout  used  was  suitable.  A  cement  slurry  with 
' Entraplast  A'  additive  has  proved  to  be  a  successful  and 
workable  grouting  material  for  sealing  well  points  into  their 
respective  locations.  If  properly  installed  a  set  cement  is 
not  required  for  continuation  of  placement  of  sand  and  seals. 

The  present  piezometer  system  although  workable  is 
not  practical  in  its  present  form.  A  more  accurate  method 
for  situating  or  seating  the  piezometer  at  the  well  point 
is  required.  The  inflatable  packers  give  rise  to  a  larger 
than  expected  volume  factor.  To  make  the  system  feasible 
the  response  time  should  be  reduced.  This  could  be  brought 
about  by  designing  a  new  packer  system  with  a  much  lower 
compliance.  Well  points  in  an  incompressible  and  imper¬ 
meable  material  will  often  result  in  long  measurement  time 
lags.  This  difficulty  could  be  overcome  with  a  system  of 
irrecoverable  pressure  transducers  when  these  devices 
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FIG.  8.4  PIEZOMETER  ARRANGEMENT  USING  IRRECOVERABLE  TRANSDUCERS 
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become  reliable,  accurate  and  economical. 

In  the  past,  coal  layers  in  the  Upper  Cretaceous 
Edmonuon  Formation  have  been  shown  to  be  free  draining  into 
the  Forth  Saskatchewan  River  valley.  Another  example  of 
this  feature  is  shown  in  the  present  study  where  irregular 
pressure  distributions  are  revealed  due  to  these  draining 
layers . 
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APPENDIX  A 

BODY  FORCE  EQUATIONS 


FIG.  A1  ELEMENT  PRESSURES  AND  RESULTANT  FORCES 

Considering  the  pressures  on  side  1-2  of  the 
element  shown  in  Fig.  Ai  the  resultant  force  F  is  given  by 

fi2  %/<*2"xi)2  +  <y2-yi>2  312 

where  d-,  9  is  a  unit  vector  perpendicular  to  side  1-2  and 

X  X 

pointing  inward. 
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By  vector  algebra  it  may  be  shown  that 


12 


r1+r2 


(yn -y0) 1  -  -x0) j 


The  sum  of  these  three  force  vectors  is  ur . 


"r  =  f12  +  p23  +  f31 


*5  L(-Riy2+R2yi-R2y3+R3y2-R3yi+Riy3)i_  (-R;jX2+R2x:rR2x3+R3x2 

-R3x1+R1x3)j  | 


The  change  in  body  force  (R  is  in  feet  of  water) 


AB.F .  =  Ay  +  u 
'  w  r 

=  -aYw7  +  "r 

=  J5YwL(-Rly2+R2yi-R2y3+R3y2-R3yi+Riy3>i 

-(2A-R1x2+R2x1-R2x3+R3x2-R3x1+R1x3)j] 

Check  on  Program  Results 

Since  the  change  in  body  force  is  identical  to  the 
seepage  force  a  check  can  be  made  on  the  program  results. 
The  following  values  are  taken  for  a  converged  case  using 
the  weaker  variation  of  permeability  and  KQ  =  2.0/  and  are 
shown  in  Fig.  A2. 
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— -  AREA  USED  IN  COMPUTING  SEEPAGE  FORCE 


FIG.  A2  ELEMENT  CONFIGURATION  AND  PRESSURES  IN  THE 
VICINITY  OF  NODE  350 


For  node  350,  the  average  gradient  in  the  x  direction  =  1.083 

» 

in  the  y  direction  =  .0225 

1  ,  2 
—  of  the  area  of  each  triangle  which  effects  node  350  =  1200  FT 

3  .  . 

The  seepage  force  F^  =  Ai^y^  =  62 . 4x 1200x 1 . 083  =  8120  lb. 

Fy  =  Aiyyw  =  62. 4x.  0225x1200  =  1685  lb. 

The  program  gives  a  body  force  F  =  8033  lb 

X 

Fy  =  1723  lb. 
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STRESS  RELEASE  FORCES 
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COMPUTER  PROGRAM 

The  following  is  a  list  of  the  important  terms  in 
the  computer  program.  Occasionally  two  notations  have 


been  used  for  the 

same  variable. 

Seepage  Section 

c  (L) 

Initial  total  allocation  of  L  storage 

locations  for  the  seepage  section. 

NUMNP 

Number  of  nodal  points. 

NUMEL 

Number  of  elements. 

NUMMAT 

Number  of  materials  (=number  of  elements) 

NUMFSC 

Number  of  free  surface  nodes. 

RO 

Density  of  fluid. 

BETA 

Free  surface  correction  factor. 

HITE, HEAD 

Factors  to  normalize  equipotential  values 

TOL 

Tolerance  for  free  surface. 

MAX  IT 

Maximum  number  of  iterations  per  cycle. 

X(  )  ,Y(  ) 

Nodal  point  co-ordinates. 

NBC  (  ) 

Nodal  boundary  condition. 

FX(  ) 

Nodal  flow  or  pressure. 

NP  (6  ,M) 

Nodal  points,  material  number  and  angle 

of  the  principal  permeability  axis 

associated  with  element  M. 

NPFS (  ) 

Free  surface  nodes. 

- 


B3 


Stress  Section 


JNUMNP , NUMNP 

Number  of  nodal  points. 

JNUMEL , NUMNP 

Number  of  elements. 

NUMAT 

Number  of  materials. 

E(  ) 

Modulus  of  elasticity. 

PR(  ) 

Poisson's  ratio. 

CO  (  )  ,  RO  (  ) 

Density. 

SX(  )  ,X(  ) 

SY(  )  ,Y(  ) 

Nodal  point  co-ordinates. 

KODE  (  ) 

Nodal  boundary  condition. 

NNP (3,1) , NP (3,1) 

Nodal  points  associated  with  element  I 

MAT  (  ) 

Element  material. 

NCORS (I) 

Seepage  element  associated  with  stress 

element  I. 

NFSC 

Number  of  free  surface  stress  nodes. 

MFP 

Number  of  slope  face  nodes. 

MBP 

Number  of  base  of  excavation  nodes. 

NCYC 

Number  of  cycles. 

DEPG 

Depth  of  grid. 

DEPE 

Depth  of  excavation. 

XKO 

Effective  stress  ratio. 

IFSN (  ) 

Free  surface  stress  nodes. 

NFP  (  ) 

Slope  face  nodes. 

NBP  (  ) 

Base  of  excavation  nodes. 
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C 

C 

c 

JL 

C 

c 

c 

c 

c 

jC 

c 


XL 

5 

c 

c 

c 

c 


c 

71 


C 

_C 

c 

c 


THIS  IS  THE  MAIN  PROGRAM  WHICH  COUPLES  THE  FINITE  ELEMENT 
SEEPAGE  PROGRAM  WITH  THE  FINITE  ELEMENT  STRESS  PROGRAM 
AN  EFFECTIVE  STRESS  ANALYSIS  IS  USED 

_ B£JZLV_LS-IL1N  I  S  M  A  O  r  .HU R_R.fi AD  LNG  _1N_JTHE_P  ERMEAB_ILJLTJ £-S  OF,  ..A - 

PREVIOUS  RUN  IF  POSITION  OF  FREE  SURFACE  NOT  STABLE 
THESE  ARE  READ  IN  FROM  A  TAPE  FILE 

THE  SEEPAGE  SECTION  HAS  A  VARIABLE  DIMENSIONING  SCHEME  WITH  A 
TOTAL  STORAGE  ALLOCATION  OF  CI21000) 

THE  STRESS  SECTION  HAS  F I XCD  DIMENSION  STATEMENTS  WITH  DIMENSIONS 
FOR  APPROXIMATELY  500  N Q DES  AND  1000  ELEMENTS _ 

COMMCN/CONTRL/HED ( 1  a )  , NUMNP, NUMEL.NUMMAT .NUMFSC . NDMP , ntype.ro, 

1  BETA, HITE, HE AD, TOL,  ma XI T , MAXMSH, NFLCO.PI  ,CU(  500)  ,CV(500) 

COM  MCN/ COMING/  E (3  )  »PP ( 3 ) ,C0(3 )  , SX ( 500 )  , SY{ 500 )  ,U(500),V<500), 
lTH(lOOO)  »STIF(1000»70) ,AP{1000) , EST I F ( 6 .6) ,ECM(3,3),EBM(3,6), 

1  ESMI3.6)  .  W  T  .  J  NU  MNP.  Jnum  EL,  J  MU  MAX* KCLQ  E  LSSlQ. )  ,  NNEX.3  x  LQ_<1QX*. - 

I  MAT (  1  000 )  , MBAND •  NEO,  M ,  LM<6),  I MX MSH ,  I  MESH ( 500 ) ,  NCORS(IOOO), 

1  I FSN ( 20  )  .NFP ( 20 ) ,  NFSC , MFP , APS (  1 00 0  )  , I H AR ,  I NCOR ( 500 ) 

1  ,XKO,NCYC  ,MBF  ,DEPG,OEPE,  NBP(20  )  ..  _  _  _ 

COMMON  0(21000) 

COM  MON /ET  S/S  I GP (  1000,4) 


READ (6*5)  I H AR 
FORMAT ( 15) 

CALL  CPLG(MAXBAN) 

C  ALL  READIN 


N1=N0+NUMMAT 
N2=N1 +NUMMAT 

_ N3-N2-*- NUMFSC _ 

N4  —  N3+NUMF SC 
N5=N4+NUMNP 

N6  =  N5  +  NUMNP  _  _ _ _ =. _ _ _ _ _ _ _ — 

N7=N6+NUMNP 

N8=N7+NUMNP 

N?=Ngt  NUMNP _ _ _ _ 

N1 0-N9+NUMEL 
N 1  l  =  M  0  +  NU  MEL *5 

_ N  1  2=N  1  1  +  NUMNP  .  _  _ _ _  _ _ _ _ _ _ _ _ 

CALL  F  CRM ( MAX 8 AN  *  C( NO)  ,C(N1  )»C  (N2 ) • C ( N  3 )  *C(N4)  »C( N5),C(N6)  , 
1_  C(  N7  )  t  C  (No  )  ,  C(  N9  )  ,  C(N10),  C(N11),  C(NI2)) _ 

CALL  CHANGE ( C( N2 ) , C ( N5 ) , C ( N6 ) . C ( N l 0 ) ) 


CALL  ASTIF 
CALL  BAND1 
CALL  STRESS 


GO  TO  71 
END 


.  MEMO  °Y  PcmxLRF-ME-NXS — 00  0430 — B-YXE-S 


» 


' 
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SUBROUTINE  CPLG(MAXBAN) 

CCNMON/CONTRL/HED ( 13), NU MNP , NUMEL, NUMMAT , NUMFSC  ♦  NDMP , NTYPE, RO» 
1BETA, HITE, HEAD, TOL, MAXI  T  ,  MAX  MSH  ,  NFL  CD  »P  I  ,  CU ( 500 )  , CV (5  00 ) 


MMAX=2l 000 
C 

c  _  _ _ _  _ _ _ 

C  READ  INITIAL  DATA 

30  REAC(5, 1006)  HED 

- 33 — R£AD(_5_,1_Q_Q0_J — NU MJ^E^NUMiLL^KUMjAT  ,  NUMFSC . NTYPE . NELCD . RO . H I TE . HEAD 

X  ,  BE  T  A  »  M  A  X  IT  ,T0L 
WRITE (6, 1026 ) 

_  WpI TE( 6 , 2000 )  HED.NUMNP, NUMEL,  NUMFSC,  RO, H I TE , HE AD , 

*  BETA , MAX  I T ,TCL 

IF  (NTYPE. EG. 0)  W R I TE ( 6 . 2 00 i ) 

_ LF  (NTYPE,  EQ,_  I)  WR  I  TE ( 6 . 2002  ) _ 

C 

C**4* *  SET  UP  VARIABLE  DIMENSIONING  ADDRESSES 

C . 

N0=  1 

Nl=NO+NUMMAT 

_ N2=N  1  F  NUMMAT _ _ _ 

N3=N2+ NUMFSC 
N4=N3*NUMFSC 

. . _  NS=  N4  +  NU  MNP  __  _ _ _  _ ■_ _ 

N6=N5+NUMNP 

N7=N6+NUMNP 

_ N8=N7-»NUMNP _ 

N9=N8+NUMNP 
N1  0— N9ENUMEL 

_  N1 1=N10+NUMEL*5  _ 

N  12  =  N1 1 +NUMNP 

c  • 

c ***** _ INPUT  MESH _ 

C 

CALL  MESH  INI MAXEAN.C ( NO )  , C { N 1  ) , C ( N2 ) , C ( N3  )  , C ( N4 )  • C ( N5 ) , C ( N6 ) • 

.X  Ct  N7  )  ,C(  NS  )  ,C(  N9  )  ,C  (Nl  0  )  )  _ _  _ _ _ _ _ _ 

IF( NDMP.NE.O )  GO  TO  30 
C 

C****  CHFCK.  ON  SIZE  OF  REQUIRED  STO  RAGE  FOR  PROBLEM  AND  COMPARE  TO  AVAIL 
C 

Nl 3  =  N1 2  +  MA  XBAN*NUMNP 

_  WPI TE (ft ,3000 )  N13.MMAX  _ 

IF(N13.GE.MMAX)  CALL  EXIT 

c 

C***»  FORMATS _ _ 

C 

1  026  FORMAT ( 1H 1  ,  *********************  SEEPAGE  INPUT  ******* ********** 

_  1  ***•//)  _  . .  . . .  _  _  .  _  .  „  . . 

1000  FORMAT ( ft  I E , 4F 1 0. 3 , I  5, F5 .2 ) 

1 006  FORMAT { 18A4 ) 

2000  FORMAT  { _ 1  P  A4  / _ _ 

1  3OH0  NUMOER  OF  NODAL  POINTS _  13/ 

2  30H0  NUMBER  OF  ELEMENTS _ 13/ 

*  30H  0  NUMBER  CF  FREE  SURFACE  CARDS,  13/  _ 


T-5TT 


. 
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A 

A 

5 


30H0 
30H0 
3  OHO 


UNIT  WEIGHT  OF  FLUID _ 

REFERENCE  FOR  POTENTIALS. 
AVAILABLE  HEAD 


E12.A/ 
El  2  .A/ 
El  2.4/ 


_ ZLQELQ  CORRECT  I  CN  FACTOR _ _ „ _ _ _ F  1  Q  .  5/ _ 

*  30H0  MAXIMUM  NUMBER  OF  I  TER  A T I  ON S . I  3/ 

7  3  0  HO  ERROR  TOLERANCE _ F10.5/1 

2001  FORMAT  ( 2  7H0  AXISYMVETRIC  FLOW  PROBLEM/) 

2002  FORMAT  (20H0  PLANE  FLOW  PROBLEM/) 

3000  FORMAT! 19H  REQUIRED  STCRAGE  =  110.  21H  ALLOCATED 

X  IIP) 

RETURN 

END 


STORAGE  = 


_  MEMORY  REQUIREMENTS  000710  BYTES 


TIT 


u  u 
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SUBROUTINE  MESHIN(MAX6AN,XK1,XK2. NPFS . AL PH A , MESH , X . Y »  NBC  * FX .ANG  , 
X  NP  ) 


SFFPAGE  M  ESEi _ CON  ELI  G  U  RAJ  1UN 


COMMON/CONTRL/HECI  18)  .  NU  MNP  •  NUMEL,  NIJMMAT  ,  NUMFSC  . NDMP , NTYPE . RO» 
1  BETA. HITE. HEAD, TOL.MAXIT, MAX  MSH,NFLCD,PI,CU(500),CV(500) 
DIMENSION  X  K  1  (  1  )  .  XK2 (  1  ) *  X (  1 )  , Y (  1)»N  BC (  1)  »FX(  1  )  . ANG ( 1  )  •  NP ( 5 • 1 ) 
DIMENSION  MESH (  1  )  , ALPHA!  1  )  , NPFS (  l  ) 


NDMP=Q 


C  *  *  *  *  READ  AND  PRINT  NODAL  INFORMATION 
C 

M=0 

_ MMM=0 _ 

60  READI5.1002)  N ,  l N , N BC < N )  ,  X ( N )  , Y ( N ) , FX ( N > 
I  F(  N.  LE  .  M  )  GC  TO  64 

MMM=MMM+ i 

MESH( MMM)=N 

IFIN.LE. 1 )  GO  TO  161 

_ NMM - N-M _ 

X  NMM  =  NMM 

DX-(X(N)-X{M) )/XNMM 

0Y= ( Y{ N)-Y(  M  )  l/XNMM  _  .  _ 

DF= (FX(N)-FX(M) l/XNMM 
MPl=M+l 

_ IF(MPl.GE.N)  GO  TO  161 _ 

NMl=N-l 

DO  61  NN=MP 1  »  NM 1 

X(  NN)=X(  NN-  1  )  +DX  _ _  _ 

Y(NN)=Y(NN-1  )  +DY 
NBC { NN ) =0 

_ FX ( NN 1=0.0 _ 

IF  (  I  M. EG. 0 )  GO  TO  61 
NBC { NN ) =NBC ( y  ) 

FX(NN)=FX(NN-1  )  +  OF  . . . . 

61  CONTINUE 

161  IF  (N.GE.NUMNP)  GO  TO  63 

_ M  =  N _ _ _ 

I  M-  I  N 
GO  TO  60 

_ 64  WPI TE(6,20  12  )  N _ 

NDMP=1 
GO  TO  60 

63  CA LL  ttRMESHl X.Y.NBC.FX) _ 

MAX  MSH=MMM 

c 

C**»*  READ  ELEMENT  CAROS  ..  _ _ _  _ 

C 

K=0 

_ M-Q _ _ 

70  RE  AD 15*  1003)  N  •  ( NP(  I  . N )  . I  =  1  . S ) • ANG ( N ) 
NP(5«N)=N 

IFIN.LE. M)  GO  TO  74  _  _ _ 


F-T5T 


' 
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OO  440  11=1.4 

DO  440  LI =1 1 , 4 

KK  = IABS(NP( Il,N)-NP(Ll,N>) 

IF  ( KK  .GT .K )  K  =  KK 

440 

CONI  I  NIUE 

MP1-M+1 

IF  (MP1.GE.N)  GO  TO  171 

NM 1  — N— 1 

00  7 l  NN=MP| ,NMl 

DO  72  1  =  1 ,4 

• 

72 

NP(  I  «NN)=NP(  I  .  NN— 1 ) + 1 

NP{ 5 iNN ) =  NN 

A  NG (NN)=ANG(  M) 

71 

171 

CONTINUE 

IF(N.GE.NUMEL)  GO  TO  73 

M=N 

74 

GO  TO  70 

WPI TE<6. 20 l 3 >  N 

NOMP= 1 

73 

GO  TO  70 

IF( NDMP.NE.O  )  RETURN 

IN=  1 

C 

75 

IM= INF49 

IF(  IM.GT.NUMEL)  I  V=N'U  VEL 

WRITE (6.200  5)  HEU ,  (  N  ,  ( NP (I. N). 1  =  1  *5)  , ANG ( N ) , N=I N, I M ) 

IF( IM.GE.NUMFL)  GO  TO  76 

IN=I M+l 

GO  TO  75 

76 

CONTINUE 

V- 

c  **** 

SET  BANDWIDTH  FOR  PROBLEM 

C 

MAXEAN  =  K  + 1 

WRI TE (6 . 2600 )  HHO iMAXSAN 

c 

c**»* 

IF(NFLCC .LE .0 )  GO  TO  78 

READ  DISTRIBUTED  FLOW  INPUT  CARDS 

c 

WR I TE ( 6  »  260 1  ) 

DO  77  N= 1  >  NFLCD 

READ! 5 . 1 004  )  I  , J, 01 J 

SI J=SQRT( ( X ( J  )-X( I)  ) **2  +  ( Y( J )— Y{  I ) ) **2  ) 

IF ( NBC(  D.NE.l  )  FX (  I  ) =FX(  I)+0.5*QIJ*SIJ 

77 

IF ( NBC (J ) .NE.l  )  F  X ( J ) =FX( J )  +  0 .5*01 J*3I J 

WRI TE(6. 2602)  I . J  *0 I J .FX(  I).FX(J) 

CONTINUE 

C 

c**** 

c 

READ  PHREATIC  SURFACE  DESCRIPTION 

78 

IF< NUMFSC .LE  .0  )  RETURN 

WR I  TEC  6, 260  3  ) 

READ ( 5 . 1 005 )  ( NPFS (  I  )  » ALPHA (  I  >  ,  1=1  , NUMFSC  ) 

c 

WRI TE (6  *  260  4  )  ( NPFS (  I  )  .  ALPHA (  I ) .  I  =  1 . NUMFSC ) 

RETURN 

rTsT 
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FORMATS 

C 

1001  FORMAT  (l5.2F10.0l 

_ I Q 0 2  F URM AT  l  1.5^_L2_.X3 . 3F  1  0.0  1 _ 

1003  FORMAT  (6I5.F10.0) 

10C4  FORMAT  (2I5.F10.0) 

_ 1005  FORMAT  (I5.F10.0)  _ 

2005  FORMAT ( 1H1 , 18A4// 

2  65H  ELMT  I  J  K  L  MAT 

_ 1 _ ANGLE//!  15.51  1  0.F1Q.3)) _ 

2600  FORMAT {  • 1  *  «  IBA4/.32H0BAND WIDTH  FOR  SEEPAGE  PROBLEM  =.I4//) 

•  2601  FORMAT ( 1 7H  DISTRIBUTED  FL0W/.46H  I  J  QIJ  FX(I) 

*  F X  (  J  )  /  )  .  _  _  ...  _  _  _  _ 

260?  FORMAT ( 2  I  5 . 3F 1 2 .4 ) 

2603  FORMAT (29H0PHREAT IC  SURFACE  DESCR I PT I ON/ . 25H 0  NODE  CORR.  AN 

*GL£ZJ 

2604  FORMAT! I  10, F15. 4) 

2012  FO  RM  AT ( 2 1 H  NODAL  CARD  ERROR,  N=  ,  13/) 

_ 2  013  FORMAT  (23H  ELEMENT  CARO  ERROR,  N=  ,.._I3/)  _ _ 

END 


s  aoocra  bytes 


rTssT 
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SUBROUTINE  WR MES H ( X , Y , NBC . FX ) 

C 

c**#*  SUBROUTINE  PRINTS  NOOAL  POSITIONS 


COM MON/ CO NT RL/ HE D ( I8),NUMNP, NUMEL . NUMMAT ♦  NUMFSC  »  ND MP  »  NTYPE. RO» 

1  BETA. HI TE. HEAD. TOL. MAX l T. MAXMSH. NFLCD.P I »CU( 500 )  »CV(  500) 

DIMENSION  X(1  ).Y(  l).NBC(l).FX(l  )  _ _ 

IN=1 

75  I  M=  I N  +  49 

_ IF  (  I  M  «GT  .NUMNP)  IM=NUMNP _ _ _ 

WR! TE (6  »  200  4 )  HED . ( N , NBC ( N ) , X ( N ) »Y(N).FX(N).N=IN.IM) 

I F(  I  V .GE. NUMNP )  RETURN 

GO  TO  75 

2004  FORMAT < IH1  ,  1 8A4// 

_ 2 _ 5  1  H  NODE  BC _ X  QBD - Y  -QBQ - E 

1  ( 2 15. 3F l 5.4 )  ) 

END 


MEMORY  REQUIREMENTS  000294  BYTES 


f 


•  3  8  5 


* 
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SUBROUTINE  READ  IN 
C 

C  STRESS  GRID  INPUT 

_ C.  THIS  SUBROUTINE  READS  AND  PRINTS  M  ATERIA1 — D  AT  A  , _ NOD  AJ — HA  I A . — Ei_£MtINJL_Q 

C  IT  GENERATES  COORDINATES  OF  INTERMEDIATE  NODAL  POINTS  AND  CALCULATES 
C  THE  BAND  WIDTH  AND  NUMBER  CF  EQUATIONS 


COMMCN/CCMING/  F  (  3  )  ,  P  3  {  3  )  ,  R  0  (  3  )  .  X(5QQ)  «._Y  l£_0  QJ^Ui^5PJ3JL,^(3jOQ-L-«- 


1 TH (  1000)  ,  S  T  I  F (1000*70)  .AP(  1000)  ,EST I F( o . 6)  «ECM(3»3),E8M<3,6)  * 

*  1  ESM(3»6),WT,  NU  MNP •  NUMEL.  NUMAT,  KODE(500)«  NP(3«1000>» 

_ 1  MAT  (  1  000  )  •  MBAND  ,  NEQ.  M  ,  LM<6),  I  MXMSH«  I  MESH  (feOO  )  ,  NCORSUOOO), 

1  IFSN { 20  ) , NFP ( 20  )  ,  NF SC , MFP , APS (  1 00 0 )  , I HAR , I NCOR ( 500 ) 

1  *  X  K  Q  »  NC  Y  C ,  M  E  F , DEPG . DE P E , NBP ( 2 0 ) 

_ DIMENSION  SHEDI  1  6) - - - - 

WRI TE(6, 700) 

C 

C  READ  PRELIMINARY  INFORMATION  _ _ 

C 

READ (5,1000)  SHED,  NU  MNP  ,  NUMEL,  NUMAT 

_ UPITF(6,2000  )  SHED,  NUMNP,  NUMEL,  NUMAT - - 

C 

C  READ  AND  WRITE  MATERIAL  PROPERTIES 

WRITEI6.2005) 

DO  10  M=l, NUMAT 

_ READ  (5,  1010)  E(M)  fi  jM  )..,_RQ  (_MJ - 

10  WRITE (6, 2010)  M,E<M),Pfi(M),KO(M) 

C 

C  READ  AND  WRITE  NODAL  DATA  AND  GENERATE  INTERMEDIATE  NODAL  DATA  _ _ 

C 

WR I T  E ( 6 • 20 15) 

_ MMM=  1 _ _ _ — - 

L=1 

READ (5,1 02  0  )  N, KODE (N)  »X(N)  »Y(N) 

_ I  MESH ( M  MM ) =N  _ _ _ 

GO  TO  40 

20  R  EAD (5,10  20  )  N, KODE ( N ) , X{ N ) , Y <N ) 

_ mmm=mmm+ 1 _ 

I MESHI MMM ) =N 
DN  =  N-L 

_ DX  =(X(N)-X(L))/DN _ 

DY  =(Y(N)-Y(L)  )/ON 
25  L=L  + 1 

_ IF(N-L) _ 50,40,30  - - - - - 

30  X ( L )  =  X(L-1)+DX 

Y ( L )  =  Y( L-l  )  +DY 

WRUE  (P.202  0  )  L,  KODE  (L).X(L),Y(L) 

GO  TO  25 

_ AO  WRITE (6 . 2020  ) _ N.KCDE(N)  ,X(N)  ,Y(N)  _ 

IFC  NUMNP— N  )  50,60,20 

50  WRITE  (6,2025)  N 

CALL  EXIT  _ _ _ i_ _ 


"TTT 


. 
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C 

C  RE  AO  AND  WRITE  ELEMENT  OATA 

C  ' 

-5.X) _ LMXMSH  =_M  MM - - - - - 

WRI TE(6. 2030 ) 

00  70  M  =  1  ,  NU  MEL 

_ REA 0(5. 1035)  M,NP(  1  , M) ,NP( 2.M » ,NP(3,M) ,MAT( M)  .TH ( M)  . NCOR  S ( M ) 

70  WRI TE(o, 2035 )  M . NP(  l  .  M ) . NP ( 2 . M)  . NP ( 3  * M )  . MAT ( M ) .  NCORS(M) 

C 

_f. _ READ  IN  ADD  I  T  I  ON  Al _ D.AXA _ 

C 

.  READ (5.1 024 )  NF SC , MF F  .  M BP . NC YC . OEPG . DEPE . XK O 

_ RE  AO (5.  1022)  (  I  FSN (  I  )  . 1= l , NF SC )  _ 

RE  AO (5.1  02  3  )  ( NF P (  I  )  ,  I  =  1  , MFP  ) 

READ(5. 1 022  )  ( N3P (  I  )  ,  I  =  1  , MBP ) 

_ WPITE(6. 1025)  NF  SC - - - 

WRI TE { 6 . 1 026  )  (  IFSN(  I), 1=1  ,NFSC) 

WRI TE(6, 1027  )  MFP 

_  WRITE16, 1026  )  ( NFP ( I ) , I = 1 , MFP )  _ 

WR ITE( 6 . 1 028  )  MBP 

WR  ITE( 6, 1026  )  ( NBP ( I  )  .  I  =  1 . MBP ) 


WR I TE ( 6 , 1031 )  DEPE 
WRITE! 6 , 1032  )  XK  G 
WRITE(6. 1033)  NCYC 


C 

c 

c 


80 


DETERMINE  BAND  WIDTH  AND  NUMBER  OF  EQUATIONS 


L=0 

DO  80  M=l,  NUMEL 
DO  8  0  .1  =  1 .2 _ _ 

11=1+1 

DO  8  0  J=  I  I  *  3 

K=  I  ARS-:  NP(  I  .  M)  — NP(  J  t  N  )  ) 


L  =  K 


IF  (K.GT.L) 

CONT I NUE 

MBAND  =  2*( L+l >  _ 

NEQ=2*NUMNP 

WR I TE( 6 .20  40  )  MBAND . NEC 


_ WRITE(6.300Qj _ : - 

IF (MBAND.LE  .  70  .  AND  .  NEQ .LE .  1  000 )  GO  TO  90 
WRITEI6.2050) 

_ CALL  EX  IT  - - 

90  RETURN 

C 

r  FORMAT  STATEMENTS _ _ _ _ _ 

C 

700  FORMAT!  1H1«  '  $$3.  STRESS  GRID  INPUT  *  11$$$$$$$$ 

1  $$$$$$'//)  -  —  —  -  - 

1025  FORMAT (///. 34H  NO.  OF  FREE  SURFACE  STRESS  NODES=,  15/) 

1026  FORMAT! 1 01 1 0) 

10  27  FORM AT ( / // . 25H  NO .  OF SLOPE FACE  NODES  =  »  15/) _ 

1028  FORMAT (/// .1 9H  NO.  OF  BASE  NODES=.l5/) 

1029  FORMAT (//. 1 5H  DEPTH  OF  GRID=.F10.l> 

1031  FORMAT (//» 2 1 H  DEPTH  OF  EXC A V AT  I ON= * F 1 0 .  1  ) _ _ 


TTsT 


'  - 
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1032  FORMAT!//,  1 4H  STRESS  RATIO=,  F6 

1033  FORMAT (//15H  NO.  OF  CYCLES=»I3) 

1  022  FO R M  AT  (  8  I  1  0  ) 

-L023 - FCrtMAXXai  l  Q  J _ 


1  ) 


///) 


1024  FORMAT (415, 3F1 0.1 ) 

3000  FORM/'T  (  *  STRESS  READ  IN  COMPLETED* 

.1  000  FORMAT ( 1 Q  A4/  316) 

2000  FORMAT!  1 OX.l 8A4 , //// 

1  1H  .  26H  NUMBER  OF  NODAL  POINTS  =  ,16/ 

- 2 — * - 2_6id — NUMB  ER_CF _ £  Lt  MENUS _ =  ,16/  _ _ 

3  1H  .  26H  NUMBER  OF  MATERIALS  =  ,16) 

2005  FORMAT!/// , 1 H  «10X»  21H  MATERIAL  PROPERTIES  // 

-  11X,  8HMAT.  NO.  »  4  X  *  9HM0DULUS  E  ,  4X  ,  l  4H  PO  I  SS  ONS  RA  T  I  O  ,  4  X  ,  1  1  HUN  I  T 

2  GHT , // ) 

1010  FORMAT (6X  ,F 1 2 .0 , 2F6 . 0  ) 

-2-Q-LQ EQRMAJL!_L±L_,_L5 >JEL7 ..0  ,  E  3  ,  F 1  7 . 1  )  _ 

FORMAT! 1H 1 , l OX, 1 8H  NODAL  POINT  INPUT  ,/// 

1 1H  » 34H  NODE  KQDE  X  COORD  Y  COORD  // ) 

FORMAT!  15,  17,  4F12.0) 


WEI 


2015 

1  020 
2020 
2025 
2 0JJL 


FORMAT ( 14 , I6,F13.3,1F12.3) 

FORMAT  (  1  HO . 2 8H  ERROR  IN  NODAL  DATA,NODE 
_EQ_R  M  X  j 1J  H  ELEMENT  DATA  ///. 


=  .14) 


1  035 
2035 
2040 


MAT 


SEEP  ELEM  //) 


1  44H  ELEM  I  j 

FORMAT  {SI5.F5.0. 15) 

FORMAT  (  14,  416,  5X , 16) 

FORMAT  1///10X,  2  2  Ft  STRESS  BAND  WIDTH  =  ,  1 6/ 

1  10X,  22H  NUMBER  OF  EQUATIONS  =  ,16) 

— H_  JEE.C  BLEM  EXCEEDS  SPEC  I  FI  ED  LIMITS 


END 


l  memory  requirements  000E62  bytes 


.. 


»  - 
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SU0R  OUT INE  FORM(MAX0AN.XK1,XK2. NPFS  *  ALPH A , ME SH , X  *  Y . NBC . F X , ANG.NP, 

X  R.C) 

C 

- — SU2R0UT  I-NE.  FQFMS  FLOW  MAXR-LCES. _ MQQ_IF-I.ES  FOR  PRESCRIBED  PRESSURE  S  . 

C***A  SOLVES  EQUATIONS  AND  COMPUTES  ELEMENT  FLOWS 
C 

_ _ _ COMMON/CONTRL/HEO (18)  »  NUMNP . NUMEL  * NUMMAT , NU MFSC . NDMP , N  TYPE ♦ RO  * 

1  BETA, HITE, HEAD, TOL.VAXI T  ,  MAX  M5H  ,  NFLCO  , PI  ,CU(500)  ,CV(500) 

COMMON /COM ING/  E(3),PR(3),CO(3),SX(500),SYC500),U(500),V(500>. 

- LXtUl  Q  03  ST IE1LQ.0  CL.  ZQ_L,  AP_L1  QHQ  I ECML3, 2±,  E3MC3 , 6J, _ 

l  ESM { 3*6  )  *WTf  J  NU  MNP .  JNU  MEL  »  JNUMAT.  KOOt ( 500 )  . NNP ( 3 . 1 0 00 )  , 

1  MAT! 1000 ), MBAND ,  NEQ,  M,  LM(6).  IMXMSH, IMESHS500),  NCORS(IOOO), 

_ 1  IFSNI20  )  , NFP ( 20 )  •  N F SC , MFP , APS (  1 00 0  ) ,  1 H AR , I NCUR I  500  ) 

1 *XKQ»NCYC*M6F,DEPG,DEFE»N0P(2Q) 

DIMENSION  XK1  { 1  )  ,XK2 (  l  ) . NPFS I  1  ) , ALPH A I  l  )  ♦ MESHI  1  )  , X (  1  )  , Y(  1  )  , NDC I  l  ) 

_ X— «_NP_L-£_.  1),P(1).C(M  AXSAN  .  1  )  _ 

CGMMON/ETS/S IGP ( 1000,4) 

DIMENSION  ANG(1).FX(1) 

_  DIMENSION  I  CAT  (  500  )  _  ....  _ _ • _  _ _ 

C 

C  ROTATE  NOOAL  VALUES  SOTHAT  CCR R ESPQN C AN CE  IS  OBTAINEO  8ETWEEN 

_ C _ STRESS  AND  SEEPAGE  NODES  IN  THE  EVENT  THAT  THEY  DO  N  O  T _ _ 

C  CORRESPOND  THE  EXCHANGE  WORKS  PROPERLY  IF  ONE  SIDE  OF  THE 

C  ELEMENT  IS  VERTICAL 

_ C . . . .  .  _ _ 

DO  7111  IX=1 ,2 

IF  { IHAR.GT.O)  GO  TO  575 

_ DO  576  JK=1,JNUMEL _ 

KK=NCORS( JK ) 

IF(KK.EQ.O)  GO  TO  576 

_ _ Nl=NNP(  1  ,  JK  )  _ _ 

N2  =  N  N  P  I 2 , JK) 

N3=NNP( 3, JK) 

_ l  1  —  NP (  1 ,KK  ) _ 

I2=NP(2,KK) 

I  3  =  NP ( 3 ,KK ) 

_ IF (  I  I ABS( N1-N2 )  .EQ . i  ) . AND.  I  I  AOS  I  I  1- 12 ) . GT.  1  )  )  GO  TO  577  .  _ 

IF {  ( I ABSI N2-N3)  .EQ . I  ) .AND.  I  I A8S ( 12-13)  . GT . 1  ) )  GO  TO  577 
IF  I  (  I ABS( N3-N1  )  .EQ . 1  )  .AND. {  I ABS ( 13- I  1 ) . GT.  1  )  )  GO  TO  577 
GO  TO  576 

577  NNP (  1  , JK  )  =N3 

NNP( 2 . JK)=Nl 

_ NNP I  3  * JK) =N2 _ 

576  CONTINUE 

575  CONTINUE 

_ 7111  CONTINUE _ 

C 

Ct***  INITIALIZATION 

DO  29  I Z  = 1 »  J  NUMNP 
29  I C  AT { I Z ) =0 

_ PI  =  3«  141592  6/180.0 _ 

NUM I T=0 

700  DO  80  11=1 .NUMNP 

DO  80  J  J  = l »  MAXB AN _ 


■ 
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80  C(JJ.II)=0.0 

DO  701  11=1 • NUMNP 


70  l 

ft  l 

R  (  I  I  )  =  0 . 6 
r n n t  imijf 

C 

IF(  IHAR.EQ.-l  )  GO  TO  2427 

GO  TO  2428 

c 

c 

2  A27 

IF  I HAR=— 1  READ  PERMEABILITIES  FROM  TAPE  FILE 

1HAR-1 

DO  2429  1=1  .  MJMEL 

READ!  1  .3319  )  XK1  (I  ) 

XK2(  I  )  =  XK 1  (  I  )  .  _  _  _  _  _  -  . 

2429 

2  A  28 

CONTI NUE 

GO  TO  2430 

CONTTNUF 

c 

IF(  IHAR.EQ.C  )  GO  TO  711 

IFINUMIT .GT  .0  )  GO  TO  612 

c 

c 

PERMEABILITY  VALUES  ASSIGNED 

nn  70  I  -  |  .  JNIJMF 1  -  - 

c 

M=NCCRS( ! ) 

IF(M.EQ.O)  GO  TO  79 

PERMEABILITY  FUNCTION  .  ....  ....  - 

79 

XKl(M)  =  0.80/(  (-APS (I  )/14A.)+20.)**3.0 

XK2 ( M ) =XK 1 (M) 

CO  NT  I  NUF 

2430 

C 

c 

CONTINUE 

OUTPUT  PERMEABILITY  . .  _ 

c 

WRITEC6. 1026  )  I HAR 

WPI TEI6.2003  )  ( XK1  ( I  )  , 1  =  1 . NUMEL >  . 

c 

c 

GO  TO  612 

ASSIGN  HOMO  AND  ISOTROPIC  PERM  TO  FIRST.  ITERATION  OF  FIRST  CYCLE.. 

c 

71  1 

DO  712  1 D  =  1  »  MJ  MEL 

XK 1  (  ID) "I  • OF—  06 

71  2 

XK2<  I D )  =  1  • OE  —  06 

CONTINUE 

rnuT  tniiif  .... 

Q  *  <• _ 

c 

c  **** 

r 

FORM  MATRICES  FOR  SOLUTION 

CALL  ODFLCWlMAXBAN.XKl  .XK2.X.Y.ANG.NP.R.C) 

IF( NDMP.NE. 0)  CALL  wRME SH ( X « Y » NBC  .FX ) 

TP/MnMP.Nir.O)  RETURN  . —  .  -  - - 

C 

c **** 

c 

MODIFY  FOR  BOUNDARY  CONDITIONS 

c 

c**** 

CALL  MOD  I FY ( M AX E AN » NBC • F X . R » C ) 

SOLVE  EQUATIONS  _  ...  -  - 

r~TST 


1 
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c 

CALL  SYMBC  <  C , R, NUMNP , VAX8AN, MAXB AN) 

C 

C 

C  PRINT  OUTPUT  OF  FREE  SURFACE 

C 

JAC-NUM  IT  +  1 

WRI TE(6 ,2349 )  IHAR.JAC 

_ ttRITJELk,  23.51  ) _ 

OO  2353  N  Z  = l  *  NUMFSC 
NZN  =  NPFS ( NZ  ) 

_ WRI TE( 6 ,23  50  )  NZN ,  XINZN),  Y (NZN )  ♦ R (NZN  ) _ _ _ 

2353  CONTINUE 

IF (  IHAR.EQ.0  )  WRITE (6,713) 

C 

c**»*  CORRECT  NOOAL  POSITIONS  ALONG  PHREATIC  SURFACE 
C 

_ NUH  I  T =NUM  IT  +  1  _  _ _ 

IF( NUMI T. GE .MAX  I T )  GO  TO  203 
ERRCJR  =  0.0 

_ DO  .  600  I  1=1  ,  NUMFSC _ „ _ I _ 1 _ _ 

M— NPFS ( f I ) 

ALP= ALPHA!  I  I  )*P  I 

_  OX  =  BET  A*R ( M) *C0S ( ALP )  /RQ  _ _ _ 

OY=BETA*RIMH>SIN(ALP)/RO 

ERROR=ERROR+OX*OX+OY*OY 

_ X(M)-X(M) FOX _ 

800  Y ( M  )  —  Y  { M)+OY 

ERROR -SORT i ERROR ) 

_ IF(  ERROR. LE.TOL  )  GO  TO  203  _ _ _ 

C 

C****  REGENERATE  MESH 

£ 

M=MESH( 1 ) 

00  fiOO  I  =2  ,  MAXMSH 

_ N=M£SH(  I  )  _  _ _ _ _ 

XNMM-N-M 

O X — (X(N)— X(M)  l/XNMM 

_ OY= ( Y(N)— Y(M)  )/XNMM _ 

MP 1 =M+ 1 

OO  601  NN=MP1 , N 

_ X (NN)-X { NN- l  ) +DX  _ 

601  Y ( NN ) =Y ( NN- 1 ) +OY 
M  =  N 

_ 600  CONTINUE _ 

GO  TO  700 
203  CONTINUE 

_ IF ( IHAR.EQ.NCYC )  GO  TO  204  _ 

C 

C  CALCULATE  THE  SEEPAGE  BODY  FORCES 

_ C _ _ _ 

DO  72  1=1, NUMNP 
C V (  I  ) =0 . 0 

_ 72 _ CUC  1  >  =  0.0 _ _ _ _ 


'4*4 
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DO  50  I=1,NUMEL 
I A=  NP (  1,1) 

I 0  =  NP ( 2  ,  t ) 

- IC  =  M£I  3  .  I.  ) _ 

TA1 23  =  X(  I C  )  * ( Y (  I  A ) — Y (18))  +X(IB)*{Y(IC)-Y(  IA))+X(!A)*(Y(IB) 
l-Y(IC)) 

_ BU=-31 .2*  <R  (  l  A)  *Y  (  IB  >-R(  IB  )  *Y  (  I  A  >  4-R  (  I6)*Y  (  IC  )-R<  .IC)  *Y  <  IB) _ 

1  +  R(  IC) *Y (  I  A )-R <  I  A) *Y ( IC >  ) 

BV=— 31 . 2* ( ♦ 1  .  *TA 1 23-R  <  ! A ) » X(  I  0) +R( IB ) *X (  IA)-R(IB)*X(IC) 

- 1 —  +  Ri-LC )  *XT_LE  )  -R  ULCJ_*X  t  I A  )  +  R  (  I  A)  *X(  IC)  ) _ 

CU(  I  A  )  =CU  {  IA  )+-0.3333*BU 
CU(  IB)=CU(  IQ  )  +  0 .3333  *BU 

_ CU(  IC)=CU(  IC  ) +-0 .33334BU  _ _ _ 

CV(  I  A ) =C V (  I  A  )+0  .33  33*BV 
CV(I8)  =  CV(IB)  +  0  .3333+BV 
CV (  I C ) =CV (  IC  H-0 .  3333*0 V 

50  CONTINUE 

I  HAR  =  I  HAR  *•  1 

_ _ _  RETURN  . .  ..  _ _ 

204  CONTINUE 
200  CONTINUE 

C 

FIND  THE  CORRESPOND ACE  BETWEEN  STRESS  AND  SEEPAGE  NODES 

DO  19  J  J=  1  ,  JNUMEL  _ ’ _ 

KK=NCORS ( J J ) 

IFIKK.EG.O )  GO  TO  24 

K1P  =  NP(  1  iKK) _ 

K2R=NP( 2 ,kK) 

K3R=NP( 3tKK) 

GO  TO  25 _ 

K1R  =  0 
K2R  =  0 

_ K3P  =  0  _ _ _ 

25  CONTINUE 
KiS=NNP( 1 , JJ) 

_  K2  S-NNP (2e JJ)  _ 

K3S=NNP< 3 , J J  ) 

IF(  ICAT(KIS)  .EO.l  )  GO  TO  26 

_ INCORIX1S )— K 1 R _ 

26  IF(  ICATIK2S)  .EO.l  )  GC  TO  27 

INCORlK2S)=K2R 

Z z _ IF{  ICAT(K3S)  .EO.l)  GO  JCL.28. _ 

INCORI K3S)=K3R 
28  I  C  AT  (  K  1  S  )  =  l 

_ I  C.ATJ  K2  SJ_=j _ 

I  C  AT  (  K  3S  )  =  1 
19  CONTINUE 

_ MCOUNT  =  0  _  .  .  _ 

DO  1119  J J=1 . JNUMNP 
KK= INCOR ( JJ ) 

_ IF(KK.EQ.O)  GO  TO  31 _ 

PSI=0. 

IF( HEAD.NE.O .0 )  PS  I = ( P ( KK  )  +RO * ( Y ( KK ) -H I T E )  ) /HEAD 
_ RKK=R<  Ki<) _ 


.=T5T 
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GO  TO  329 
31  RKK=0. 

PSI=0.0 

- 32-9 - CONTINUE _ _ _ 

C 

C  PRINT  AND  PUNCH  OUTPUT  OF  PRESSURES  AND  POTENTIALS  AND  OTHER  DATA 

_ C _  _  _ 

MCOUNT=MCQUNT— 1 
I  Ft MCOUNT .GT  .0 )  GO  TO  119 

- MCOUNT  =  50 _ _ _ 

WRITEI6.2008)  HED 

119  WRITE(6,5)  J J ,  INCOR< JJ  )  , SX( J J  )  . SY ( JJ ) ,  U(JJ>,  V ( J J  )  , RKK . PS l 

_ _ PUNCH  2377, JJ ,  INCOR< JJ ), SX (JJ)  ,SY ( JJ)  ,U ( JJ)  ,V( JJ  ), RKK. PSI . 

1  X ( K  K )  , Y ( KK ) 

1119  CONTINUE 

_ DO  1  20  JK  =  1  .  JNUMFL _ _ _ 

KK=NCORS ( JK ) 

N1=NNP( 1 , JK) 

_  N2— NNP ( 2 , JK )  _  _ 

N3=NNP( 3 . JK) 

I  1 =NP (  1  ,KK) 

_  _ L2  HN£LL2_tHl<Ll _ _ 

I3=NP( 3 ,KK ) 

XAVE-{ SX( N1  )  +  S  X  (  N2  )  +  S  X  (  N  3  )  ) *0.3333 

_  YAVE=(SY (N1  ) +  SY( N2 ) +SY( N3) ) *0.  3333 

PUNCH  2376 . JK.Nl  ,N2, N3. KK ,  I  1  ,  12 , 13. XAVE ,  Y A V E . X K l  ( KK  )  , S I GP ( JK .  1  )  . 

X  S I GP ( JK  «  2 )  »  S IGP( JK, 3)  , SIGP ( JK. 4) , APS( JK ) 

_ 120  CO  Nil NUE _ 

c 

C  READ  PERMEABILITIES  INTO  TAPE  FILE 

_ C _ _ _  ...  _ _  _ _ 

WRITE!  1 ,3319  )  (  XK 1  (  I  )  ,  I -1  ,  NUMEL  ) 

C 

_ C  *  *  * ♦  SOLVE  FOR  ELEMENT  FLOWS _ 

C 

CALL  ELFLOW  (  XK  1  ,  XK2  ,  X,  Y,  ANC-.NP,  R  ) 

32  WR I T  E ( 6 ♦ 20 15) 

STOP 

C 

C****  OUTPUT  FORMATS 
C 

_ 1026  _ FORMAT!/  11H  CYCLE  N0.=  ,I3)  _  _ _ 

713  FORMAT ( 4SH  INITIAL  ITERATION  WITH  HOMOGENEOUS  PERMEABILITY) 

2  376  FORMAT  (8I3.2F6.1.E10.3.3F7.0.F5.0.F8.0) 

_ 2  377  FORMAT  (2I3.2F7.2.2F10.1  .2F7 . 2 .2F8.2  ) _ _ 

5  FORMAT (5X, I3.13X, I3,9X,F7.2»4X,F7.2,3X.F10. I ,3X,F10.1.7X,F7.2, 

1  8X.F7.2) 

_ 3.3  1  9  FORM  AT  (  6  1  2 .4  )  _ _  _ _ _  _ ; _ 

200 J  FORMAT  { 28H0M AT ER  IAL  PERMEABILITIES// 

1  (  1  0  E 1  1  .  3  )  ) 

_ 2349  FORMAT (  1  HI  ,22H  FREE  SURFACE  LOC AT  I  ON . 20X  ,  1  1 H  CYCLE  N0.=  ,I3. _ 

1  J5H  ITERATION  N0.=,I3/) 

2351  FORMAT ( 5X»  5H  NODE , 5 X  ,5H X-ORD . 5 X . 5H Y -ORO , 1 2 X ,  8HPRESSURE/) 

2350  FORM  AT  {  5X  »  15.  2F10.4,  E20.5)  _ _ _ _ _ _ 


535 
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2008  FORMAT (  1 H 1  .  1  8 A4// 1 2 H  STRESS  NODE.3X.13H  SEEPAGE  NODE  »  6  X . 

1  6H  X-0R0.SX.6H  Y-0KD,5X,8H  X  FORCE, 5X.8H  Y  FORCE, 5X, 

1  9H  PRESS  OR  E , 5X  »  10H  POTENTIAL) 

- 2-01.1 — EOR'AA-TXSX  .15-. 2F10.4.  ?F  ?D.S> _ _ 

2012  FORM  ATI  I  5 ,4F 1 0. 2 ,F 1 0 . 0, I  5) 

2015  FORMAT  (15H1END  CF  PROBLEM/) 

_ 2016 _ FORM  AT_  {  615,  E15.5  ) _ _ _ 

END 


-EU  jLS — O  0  20  92  BYTES 
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C 

C 

£ 

C 

C 

C 

C 

C 

JL 

C 


42 

C 

C 

C 

c 

c 


27 


60  1 


6C0 
C _ 

c 

c 

c 


SUBROUTINE  CHANGE! NPFS. X, Y,NP) 

THIS  SUBROUTINE  TAKES  THE  BODY  FORCES  CALCULATED  IN  SUBROUTINE 
FORM  AND  ASSIGNS  THEM  TO  THE  PROPEP  STRESS  NODES 

_ IT.  CHA NG FS  THE  TNTFk i\A! NODAI RD Si. OJ3NS_JTU UHL RQSJLXI ON  QE _ 

THE  FREE  SURFACE  FOUND  IN  FORM  AND  GENERATES  NEW  POSITIONS 
IT  ALSO  CALCULATES  THE  RELEASE  FORCES  FOR  THE  SLOPE  FACE  DUE  TO 

_ THE  READJUSTMENT  OF  NODAL  POSITIONS  _  _  _ 

IN  ADDITION  TO  THE  EASE  NODE  FORCES 

THE  SLOPE  FACE  NODES  ARE  READ  IN  GOING  UP  AND  THE  EASE  NODES 
_ ARF  RFAD  IN  FROM.  THE  LFFTTfl  THF  SLOPE  FACE _ 

CCMMON/CONTRL /HE D (  1 8)  .NUMNP, NUME  L.NUMMAT . NU MFSC « NDMP . NT Y PE • RO, 
IBETA.HiTE.HEAO.  TQL.PAXI  T.MAXMSH.NFLCD.PI  ,CU(  500)  ,CV(500) 

COM MON/ COM  I NG/  E(3),FR(3),CO<3).SX(S00).SY(50O),U(500),V(500), 
1TH{  10  00)  «STIF(  1000*70)  .AP(1000)  . EST I F ( 6 « 6 ) . E CM { 3 , 3 )  *EBM(3.6). 

1  ESM I  3  *  6  )  *  WT  *  J  NU  MNP  «  JNUMEL  «  JNUMAT  *  _ 

1  MAT(  1000)  ,MBAND,  NEQ  *  M,  LM(6),  I MX MSH ,  I  MESH ( 50 0 ) ,  NCURS(IOOO), 
1  IFSN(20  ) »NFP(20 ) •  NF SC * MFP *  APS (  1  00 0  )  •  I HAR *  I NCOR { 500 ) 

1  ,XKO,NCYC,M3F.DEPG,DEPE.NBP(20)  .  _ _  _ _ _ _ 

DIMENSION  I C  A T ( 5 00  ) 

DI MENS  I  ON  X (  1  )  , Ytl  ) »  NP (5*1  ) • NPF  S ( 1  ) 


DO  42  I  =  1  »  J  NUMNP 
U(  I ) =0 

V(I)=0  _  _  _ _ _ 

CHANGE  THE  POSITION  OF  THE  FREE  SURFACE  STRESS  NODES  AND 

REGENERATE  NEW  INTERMEDIATE  LOCATIONS _ 

MODIFY  FOR  BOUNDARY  CONDITIONS  WHERE  REQUIRED 

DO  27  1  =  1. NU  MF  SC  _ 

NA=NPFS {  I  ) 

NG= I FSN (  I  ) 

SX { NB )  =  X( NA  ) _ 

S Y ( NB ) = Y ( N A ) 

CONTINUE 

_ M= I  MESH ( !  )  _ _ 

DO  600  1  =  2 .  I  MX  MSH 

N=I MESH(  I  ) 

X  NM  M  =  N—  M _ 

DX=(SX(N)-SX(M)  ) /X  NM  M 
DY=(SY(N)-SY(M)  )/XNMM 

_ MR  1  =  M  +  1  _ _ _ _ _ 

DO  601  NN=MP1 , N 
SX(NN)  =  SX(NN-1  ) +DX 

SY!NN)  =  SY(N N - 1  ) +CY _ 

M  =  N 

CONTINUE 

CALCULATE  RELEASE  FORCES  FOR  HORIZONTAL  STRESSES 
MODIFY  FOR  BOUNDARY  CONDITIONS  WHERE  REQUIRED 


DO  39  1=1. MFP 

IF (  I  •  EO  .  1  )  GO  TO  20 

IF  (  I  .EG .  MFP)  GO  TO  30 
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GO  TO  40 
20  M1=NFP( I ) 

M2=NFP( I  +  l  ) 

- DM=DEPG-SY ( Ml  ) _ 

DMP 1 =DE  PG  — S Y ( M2 ) 

U(Ml)=-;:Q(3)*XK.a*0.125*<3.*<DM**2)-2.*DM*DMPl-DMPl**2) 

_ IF(KC0E(M1 ) .EQ.l 1 )  U(MI)=0.0 _ _ 

GO  TO  39 

30  Ml=NFP(I-l) 

_ M2-ENFP  (  I  ) _ 

DMM  l  =  DEPG-SY( Ml  ) 

1  U{  M2  )- —  C  0(  :*  )*XK0*0.125*(  OMM  1**2) 

_ 1F(K00E(M2)  .EQ.l  1  )  U(M2)=0.0  _ _ _ _ _ _ _ 

GO  TO  39 

40  M 1 =N  FP {  I — l  ) 

_ M2  =  NFP(  I  ) _ _ _ 

M3=NFP ( I+l ) 

OMMl=OEPG-SY(Ml ) 

_ DM=DEPG— SY  (  M2  )  _ _  _  _ _ 

OMP1=OEPG— SY(M3) 

U( M2) =— C0(  3  )*XK0*0.125*{ OMM 1**2  +  2.*DMM1 *DM— 2 • *D  M*D  MP 1  — O  MP l **2 ) 

_ _ IF(K0PE(M2)  .to.  I  1  )  U  (  M2  )  =  0  »  0 _ __ 

39  CONTINUE 

C 

_ C _  CALCULATE  RELEASE  FORCES  FOR  VERTICAL  STRESSES  _ 

C  MODIFY  FOR  BOUNDARY  CONDITIONS  WHERE  REQUIRED 

C 

_ DQ...A3  1  =  1.  M£1E _ 

I  F(  I  .EQ  « 1  )  GO  TO  44 
IF (  I  *  EQ  *  MBP)  GO  TO  45 

_ GO  TO  46  _ _ _ 

44  M l =  N6P ( I ) 

M2  =NBP (  I +1  ) 

_ V(  Ml  ?=  0(  ;  )»DEPE*(  SX(M2)-SX(M1  )  )  -+  Q  .  -3 _ 

I F ( K  00  E(Ml).EQ.ll)  V(M1>=0.0 
GO  TO  43 

_ 45  _ Ml  =NBP  (  I  )  _ _  _ _ 

M2  =  NBP(  1-1  ) 

V<  Ml  )  =  CO  <  3 )*DEPE*( SX( Ml  )-SX{ M2 )  )  .5 

_ IF ( KODE ( m 1 )  .EQ  .11)  VIM  1  )  =  0. 0 _ 

GO  TO  43 

46  M  1  =  Nf3P  (  I  —  1  ) 

_ M2=NBP(I)  _ 

M3=NB°( I+l ) 

V(  M2  )  =  C  Q  (  i  i*DEPE*(  SX(  M3  )-SX  (  Ml  ))*  C  .  5 

_ IF ( KCOE ( M2 )  .EQ . 1  1 )  V(M2)=0.0 _ 

43  CONTINUE 

WR  I TE ( 6  » 1 30 )  I HAR 

_ WR ITE( 6 i 47  )  _ 

DO  53  1=1. ME  P 
NI =NFP(  I  ) 

_ 53 _ WR  I  TE  (6.48  )  NFP(I).U(NI) _ 

WRI TEI6.49) 

DO  54  1=1. MBP 

M I =  NBP (  I  )  _ 
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54  WR I  T  E ( 6  »  50  )  NBP(I),  V(M!) 

C 

C  ASSIGN  BODY  FORCES  FOUND  TO  THE  PROPER  NODES 

_ C _ MODIFY  FOR  BOUNDARY  CONDITIONS  WHFRE.  REQUIRED _ : _ 

C 

DO  36  I  =  1 • J  NUMNP 

_ 36 _  I  CAT  (  I  )  =0  _ _ 

DO  37  I  =1  ,  J NUMEL 
M=NCCRS( I ) 

_ IF  (  M.EQ.Q)  GO  TG  3  7 - -  - 

.  I  A  =NNP {  1,1) 

I  Ft=NNP  (2.1) 

_ _ I  C=NNP  {  3,1)  _ _ _ 

[MA=NP( 1 ,M ) 

I  MB  =  NP ( 2  »  M ) 

_ LM  C^JSLP  (  3  .  M  ) _ - — - 

IF (  IC AT (  I  A )  .EQ . 1  )  GO  TO  32 
IF(KOOE( IA ) .EO.I 1 )  GC  TO  32 

_ IF(KOOE(  IA) .EQ. 10)  CU(IMA)  =  0.0  _ _ _ 

IF(KODE(  IA )  .EQ.01  )  CV(IMA)=0.0 

Ul  I  A  )  =  U (  I  A ) +CU(  IMA ) 

_ YJ_LA_Lf  Ml  J.  AJ_i  CV.LIM  AJ _ : _ 

32  IF  (  ICAT(  IB)  .EQ.  1  )  GC  TO  33 
IF(KOQE( IB) .EG. 1 1 )  GC  TO  33 

_ _  IF ( KODE <  I  8)  .EG . 1 0 )  CU(IMB)=0.0  _ 

IF{ KODE( IB) .EG. 01 1  CV(IMB)-0.0 

U(  IB) -U<  IBMCU<  IMS  ) 

_  V(  IP)-V(  I  9  )  -*■  C  V  (  I  MB) _ 

33  IF(  ICAT(  IC)  .60.1  )  GO  TC  34 
IF (KOOE(  IC)  .60. 1  1  )  GO  TO  34 

I F ( KODEi  IC )  .EQ  .  1 0 )  CU(IMC)=0.0  .  _  _ 

IF(KODE(  IC)  .EG. 01  )  CV(IMC)=0.0 

U<  I C )=U{  IC ) +CU (  I  MC ) 

_ v  1 1  c )  -  v  1 1  c  L±.CV..(.  I  vc J _ 

34  I C  A  T (  IA )  =  1 
I  CAT (  IB)  =  1 

_ I  C  AT  (  10  =  1  _ _ 

37  CONTINUE 

RETURN 

_ C _ _  ....... - * - 

C  FORMAT  STATEMENTS 

C 

130  FORM  AT (  1H1  ,  1  1 H  CYCLE  N0.=  ,  13/)  _ 

47  FO  RM  AT (  18H  SLOPE  FACE  FORCES/ 

l  5X,  5H  NOD  E  » 1 0  X  »  3H  X  FORCE/) 

_ 4J3 _ FORMAT  t  7X  .  I3.10X.F10.1  ) _ 

49  FORMAT ( // »  17H  EASE  NCOE  FORCES/ 

1  5X  *  SH  N0DE.10X.eH  Y  FORCE/) 

_ 50  FO  R  M AT (7X. I 3» 1 0X.F1 0.1  )  _ _ _ 

ENO 
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SUBROUTINE  QDFLOW ( MAXEAN » XK 1 , XK2 , X , Y , ANG,NP,R,C> 

C  \ 

C***4  SUBROUTINE  COMPUTES  FLOW  MATRICES  FOR  QUADRILATERAL  OR  TRIANGULAR 
_ C*  »  »  * _ FI  ■  FMFNTS  AND  Anns  Tfl  GIORAI  F1QW  MATRICES _ 1 _ 

c 

common/ccntrl/hed ( i a> . numnp, numel.nummat ,numfsc . ndmp, ntype. ro. 

_ 1  BE  TA  ,  HITE.  HEAD.  TOL,  MAX  I T ,  MAX  MSH  ,  NFLCD.P  I  , CU ( 500  >  . CV (  5 00 >  _  ... 

COMMON/ EL MT /  XX(5).YY{5)  . S (5  *5)  ,GG(5),XKll » X  K22 ,  XK1  2 . NN 
DIMENSION  XK 1 ( 1 ) . XK2( 1 ) , X ( 1 ) , Y ( 1 ) , ANG( 1 ) ,NP( 5. 1 ) ,R ( 1 ) »C( MAXE3AN, 1 ) 

_ DO  3Q0  NN  =  1  ,  Ml  MFL _ 

C 

C***«  INITIALIZE  MATRICES 

DO  20 l  11=1,5 
GG(  I  I  )=0 .0 

_ DO  201  J J=1 . 5 _ 

201  S(I I . JJ)=0.0 

XK3=ANG( NN)*PI 

_ _ HAT  =  NP(  5  (  NN  )  _ _ _ _ _ _ 

CC=COS (  XK3  ) 

SS=S  INI XK3 ) 

_ XK  1  I  =  XK  1  I  MA T  )  * CC* C CJLXK2XMAT)  *SS*SS _ 

XK22=XK2( MAT ) *CC*CC+XK1 ( MAT ) *SS*SS 
XK12=SSACC*(XK1 t  MAT )  — XK2 (MAT)  ) 

XX { 5  >  =0 . 0 
YY(S)=0.0 

_ DO  100  J  =  l,4 _ 

I =NP ( J , NN ) 

XX ( J ) =X {  I  ) 

YY<  J)  =  Y(  I  )  _ _ 

X X ( 5 ) -XX( 5  )+X(  I  )/D 
100  YY(5)=YY< 5)+Y( I )/0 

_ IFINPI3.NN)  .EC.NP(A.NN)  )  GO  TO  110 _ 

C 

C ****  FORM  QUADRILATERAL  FLOW  MATRICES  FROM  TRIANGLES 
MM- 4 

CALL  TR  I FL (  1,2,5) 

_ CALL  T  E  I  FL (2,3,5) _ 

CALL  TRIFL(3,4,5) 

CALL  TRI FL (4,1,5) 

_  IF( NDKP.NE.O)  GO  TO  300  _ 

C 

c *  +  **  reduce  center  node 

_ Q _ — - 

DO  200  11=1,4 
COM  =  S (  I  I  ,  5  )  /S ( 5  »  5 ) 

GG(  I  I  )  =GG ( J  I  )  — C  C  M*GG ( 5 ) _ 

DO  200  J  J— 1 »  4 

20C  S(II,JJ)=S(II  , JJ )-CCM*S(5  ,  JJ  ) 

_ GO  TO  250 _ 

110  CONTINUE 
C 

c****  FORM  MATRICES  FOR  SINGLE  TRIANGLE  _ . 
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C 

MM=3 

CALL  TR IFL( 1 , 2. 3) 

IF( NDMP.NF.O)  GO  TO  300 

c 

c*  *** 

c 

ADD  ELEMT  MATRICES  TO  GLOBAL  MATRICES 

250 

CONTINUE 

DO  202  11=1. MM 

J_1  =  NP(  II.NN1 

R (LI ) =H( LI ) 4GG( I I > 

DO  202  J J=1 , MM 

K2=NP< J J , NN )-L 1 +1 

20  2 

IFIK2.LE.0)  GO  TO  202 

C ( K2  .  L 1  )=C(K2,L1  )4S(I  I.JJ) 

CONTINUE 

300 

CONT  INUE 

RETURN 

END  .. 

_  MEMORY  REOU I REMENTS  0006B0  BYTES 
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C 

c*  *** 

c 


SLDROUT I NE  TRIFLI  I • J«  K ) 

SUBROUTINE  COMPUTES  FLOW  MATRICES  FOR  TRIANGULAR  ELEMENT 


COMMON/CONTRL/HED { 1  a )  , NUMNP, NUMEL.NUMMAT ,NUMFSC , NDMP . NTYPE. RO. 

1  BETA, HITE, HEAD, TOL, MAXI  T , MAX MSH , NFL CO , P I  , CU( 500)  ,CV( 500) 

_ CCMMON/ELMT /  X X ( 5  )  , Y Y ( 5 )  ,  S ( 5 , 5 )  , GG ( 5 ) *  X K 1  1 , X K22 * XK 1 2 , NN  _ 

A 1  =  0 . 5 

P21 =YY ( J )-YY( K ) 

P2?=YY ( K  l-YYl  I  ) 

P23=Y Y i  I  ) — Y  Y ( J ) 

P3  l--XX(K)— XX(  J) 

P32  =  XX (  I  ) — X  X ( K )  _ _ 

P33  =  X X ( J )  —  XX (  I  ) 

D=P33*P22- P2 3* P3 2 

_ I  F  (  NO  MP.NE  .0  )  RET  UR  N _ 


IF(D.LE.O.O)  GO  TO  500 

IF ( NTYPE. ME.  1)  Al  =  (  XX( I  )  +  XX( J)  +  XX(K) )/6.0 

T1 1=XK  1  1  +  P2l+XK12*P31  .  _ _  _ 

T12=XK1I*P22+XK12*F32 
T 1 3=  XXI 1*P23VXK12  *P33 

-T2  1  =  XK,  1 2 P 2. 1  •*- X < 2 ?  »P3  I _ 


T22=XK1 2*P22+XK22*F32 
T23=XK12*P23+XK22*P33 
COM=A i/D 
S( I  .  I  )  =  S (  I  ,  I 
S  (  I  • J  >  =S (  J  ,  J 
S(  I  ,K)-S1  I,K)+CCM*(P2l*T13+P3  i  *T23  ) 


S(J,J)=S(J,J)+CCM*(P22*T12+P32*T22) 
S(J»K)=S,J»K)  +C0V$(P22*T I  3 +P32*  r 23) 

S(  K.K)=S{K,K)  +  CCM*(  P23*T1  3*P33*T23)_. 
S  (  J  »  ! )=S(  I  ,  J  > 

S  (  K  ,  I  )  =  S  (  I  ,  K  ) 

S(K,J3=Si J.K) _ 


500 


C0M=A1 *RO 

GG{  I  ) =GG (  I  )  — T2I *CCM 
GG(J)=GG(J)-T22*CCM 
GG ( K )  —  GG ( K  )  —  T23*C0M 
RETURN 

WRITE! 6 , 3000  )NN 


NDMP  = 1 
RETURN 


C 

C**** 

c 

3000 


ERROR  MESSAGE 

_FOR  MAT  (  3  OH  ZERO  OR  NEGATIVE  AREA  ELEMT,  13/1 


END 


L  MEMORY  . REQUIREMENTS  0006A4  BYTES . . 
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SUBROUTINE  MOD  I FY ( M M , NQC . FX . R. C ) 

C 

+  *  SUBROUTINE  MOCIFIES  EQUATIONS  OF  NODES  WITH  PRESCRIBED  PRESSURES 

- OR  FI  riw S _ _ _ 

C 

COMMON /CONTRL/HEC (13), NUMNP , NUMEL , NUMMAT , NUMFSC , NDMP , NTYPE, RQ» 

_ 1  0ETA»HITE»HEAC*T0L»  WAX  I T , MAX  MSH , NFL  CD  »  P I , CU (500)  »CV(SOO) 

DIMENSION  FX ( 1 ) ,R ( 1 ) , C ( MM , 1 ) , NBC ( 1 ) 

DO  70  N= 1 , NUMNP 

- LfLLN  BC(  N)  .  G  T  .  0  )  GO  TO  72 _ 

•  R(N)=R(N)+FX(N) 

GO  TO  70 

_ 72  N1  =  N-H  ______  _  _  _ 

DO  73  K=2.MM 
N2=N1 -K 

_ LE _ t.N Z.LE.O)  GO  TO  74 _ 

R(N2)=R(N2)-C(K,N2)*FX(N> 

C ( K  »  N2 ) =0 • 0 

_ 74  V  =  N  +  K-l  _ _ _ _ 

IF  (M.GT. NUMNP)  GO  TO  73 
R(M)=R(M)-C(K,N)*FX(N) 

I _ C  (  K  t  N  )  -  0  .  O _ _ _ _ _ 

73  CONTINUE 
C{ 1  ,N)  =  0.0 

_  R(N)=FX(N)  _  _  _ _ 

70  CONTINUE 
RETURN 
END 

_  MEMORY  REQUIREMENTS  000376  BYTES 
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SUBROUTINE  S Y MBC I  A , B . NK , MB , MM AX  ) 

C 

C***o  SOLUTICN  OF  SYMMETRIC  EANOEO  EQUATIONS  IN  SINGLE  SUBSCRIPT  ARITH. 

- CJtJkJt  A — SY  MMEIRLC.  . BAND  N?  ATJ1 1  £L  JW  TlhLJlI  MEN  ST  QNS iU.MMAXiTiN-I.-lS  _  SIN  RED  IN _ 

C*  ♦  *  *  VECTOR  A. 

C 

_ DIMENSION  A (  1  )»R ( 1  ) _ 

MB1 =  MB— 1 
NNN=NN— 1 

C 

C ****  TRI  ANC-ULARICE  MATRIX  BY  GAUSS  ELIMINATION  WITHOUT  PIVOTS 
C 

_ I  1=1  _ _ ; _ i _ 

DO  300  N  = l  ,NNN 
CC=  A (  I  I  ) 

_ L£X£.CUL£LlQ_«  Q  1.  GQ  ..TQ..250 _ 

Jl= I  1  +  1 
J2= I  I +MB 1 

_ NE=NN-K  _ _ _ _ _ 

IFINE.LT. MB1 )  J2=II+NE 
M=I 1-1 

_ Qg_200  J  =  J1  » J2 _ * _ 

M=M+MMAX 
C= A ( J )/CC 

_  IF(C.EQ.O.O)  GO  TO  200 _ 

K=M 

DO  100  I  =  J  »  J  2 
K  =K+  1 

100  AIK )  =  A ( K ) — C  * A (  I  ) 

A  I  J  )  =  C 

_ 2  00  CONTINUE _ _ 

250  CONTINUE 

I  1= I I+MMAX 
300  CONTINUE 

C 

C  *  *  *  *  REDUCE  FORCE  VECTOR 

_ C _ 

1 1=  1 

DO  500  N=1  .NMv 

_ CC=  A (  I  I  ) _ 

IF( CC.EQ  .0 .0  )  GO  TO  450 

J1=I 1+1 

_ J2= I  I +MS1 _ 

NE=NN—  N 

IF( NE.LT.MB1 )  J2=II+NE 

_ C  =  8  < -Nl _ 

L  =  N 

DO  400  J=J1 » J2 

_ _ _ L=L+1  _ _ _ _ 

400  B(L ) =B(L )— A( J  )*C 
BIN) =C/CC 

_ 450  CONTINUE _ 

I  1= I  I +MMAX 
500  CONTINUE 

CC- A (II)  _ 
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IFCCC.NE.0.0  )  0(NN)=6{NN) /CC 
C 

C****  BACKSUBST ITUTE  ECUATIONX 

c 

N=NN 

I I=MMAX*CNN-2 )+l 

_ DO  700  I =2 ♦ NN _ 

N=N-1 

I  F( AC  1 1  )  .EQ  .0 .0 )  GO  TO  650 

_ J  1  =  I  I  +  1 - - 

J2= I  I +M81 
NE=NN-N 

_ I F ( NE.LT.MB1 )  .  J2=  I  l  +  NE _ 

C=B ( N ) 

L  =  N 

_ nn  600  J  =  J1  . J2  - 

L=L+l 

600  C=C— AC J)*Q(L) 

650  CONTINUE 

!  1=  I  l -MMA X 

_ 700  CONTI  Mile - - - 

RETURN 

END 


L  MEMORY  REQUIREMENTS  00062C  BYTES 
48.22  RC=0 


PAGE  0002 


-'•VST" 
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SUE? ROUT  INE  ELFL0W(XK1,XK2»X,Y,ANG,NP,R> 

C*++*  SUBROUTINE  COMPUTES  FLO*  VELOCITIES  IN  ELEMENTS 
C 

_ CQMMON/CQNTRL/HED  (  18).  NUVNP.  NUMEL.  NUMMAT  .  NU  MF.SC_jJSD.MP  t,N.T,Y_EE-t  R-£b - 

1  BETA, MITE, HE AD, TOL. MAXI  T , M aX MSH , NFL CO ♦ P  I  ,CU(500)  ,CV(500) 

DIMENSION  XK1I1),  XK2 (1).X(1),Y(1),ANG(1).NP(5.1).R(1J,YY(6) 

_ MCOUNT  =  0  _ _ _ _ _ 

DO  400  N  =  l  , NUMEL 

IF(NP(3,N) ,EQ.NP(4.N)  )  GOTO  200 

X  M  -  O  .  .  

YM=0  . 

.  RM=0.0 

_ DO  100  J=l  ,  4  _ _ 

I =NP ( J , N ) 

RM=RM+0 . 25*R ( I) 

_ 

100  YM- YM+0 * 25*Y ( I) 

MM=4 

_  GO  TO  250  _ _ _ _ 

200  XM=  0,0 
YM  =  0.0 

_ RM  -  0.0 _ 

DO  500  J= 1 • 3 
I  =  NP(J,N> 

XM  =  XM  +  XII)  /  3  •  0_ _ : _ 

RM  =  RM  +  R ( I )  /  3.0 

500  YM  =  YM  +  Y{ I )  /  3.0 

MV  =  3 

250  CONTINUE 

MAT-NP { 5  »N) 

XK3=ANG(N )*PI _ 

CC=COS( XK3 ) 

SN=  S I N ( XK  3 ) 

_ Y1K=CC*XM+SN*YM _ _ _ 

Y2K=CC*YM-SN*XM 
01  =  0 . 0 

_ 02=0.0 _ 

Dl  =  0 .0 
C 

C  t ***  LOOP  ON  ELEMENTS  FOR  QUADRILATERAL  ELEMENT _ 

C 

J=NP (4 ,N ) 

_ DO  30  0  NN  =  l , MM _ 

I=J 

J=NP ( NN iN) 

_ Y 1  I=CC*X(  I ) +SN»Y(  I  ) _ 

Y1J=CC*X(J)+SN*Y(J) 

Y2  I  =  CC*Y (  I  ) -SN* X {  I  ) 

_  Y2 J=CC*Y I J )-SN*X( J ) _ 

YY ( 1 )  =Y2 J-Y2K 

YY 12)  =Y2K-Y2I 

_ YY  t  3  )  =Y  2  I  —  Y  2  J _ _ _ 

YY { 4  )  =Y 1 K— Y  1  J 
YY(  5  )  =Y  1  I-Y  IK 

YY( 6)  =  Y  1  J  — Y 1  I  _ 
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D=YY< 6>*YY( 2>-YY(3)*YY(5> 

01=01+0 

01=01— XK.1  ( MAT  )*(YY(  1  > ♦ R  C I  ) ♦ YY ( 2 ) ( J ) +Y Y( 3 )*RM+SN*RO*Ol 

_ G2=Q2--X.K2  (  MA  T  )  »  (  YY  (  4  1  ♦  R  (  I  )  4  Y  Y  (S  1»R1^114YYL6 1WHM4CC »HQ*DJ - 

300  CONTINUE 
C 

_ C  *  *  **  OUTPUT  FLOWS  _ _ 

C 

01=01 /D1 

_ o?=Q2/n  i _ _ _ _ 

Q3=S0RT(Q 1**2402 **2) 

Q4=ATAN2( Q2.Q1 >/PI+ANG(N) 

_  ‘ _  MCOUNT  =  MCOUNT— 1  _ _ _ _ _ 

I F ( MCOUNT . GT .0  )  GO  TO  350 
MC  OUNT  =  55 

_ WRITE(6.2Q00)  HEP _ 

350  WRITEI6.2001  )  N  ,  X M  ,  YM , 0 1 . 02 , A NG ( N ) , 03 , 04 
400  CONTINUE 

RETURN  _  _  .  _ 


C *+**  FORMATS 


2000  FORMAT (  1H1  •  18A4//5X.5H  ELMT < 5X» 5HX—  ORD  »  5  X  »  5HY— CRD . 1 4X , 6H 1 -F LOW . 
X  14X . 6H2-FLOW, 1  OX .5HANGLE ♦ l OX . 1 OHTOTAL  FLOW . 6X .9HD I  RECTI  ON/  ) 

2001  FORMAT ( 5X. 15 .  2F10.4.  2E20.5.  F15.4.  E20.5.  F15.4  )  _ 

END 

MEMORY-  RHQLL  LREMFNTS  0  0  Q.7.B-4 — BYT  ES - 


' 


B31 


I  IV  G  COMPILER  ASTIF 


03-29-72  23:00.58  PAGE  0001 


SUB ROOT  I NE  ASTIF 
C 

C  THIS  SUBROUTINE  TAKFS  EACH  ELEMENT  IN  TURN  AND  FORMS  THE  ELEMENT  S' 

£ - MAJJLL-X LBY  CA  LL  I N  G EL  S.T  I  Fl.TT  ASS  FMRI.FS  THF  Fl.FMF  NT  STIFFNES  S  E  S  IN' 

C  KSTIF  ,  ASSEMBLES  THE  APPLIED  LOAD  VECTOR  (AP),  AND  MODIFIES  THE 

C  ASSEMBLASES  FCR  DISPLACEMENT  BOUNDARY  CONDITIONS  (BY  CALLING  MOD  I F* 

c _ _  _ 

COM  MON/C  CM  INC/  E  (  3  )  .  PP  (  3  )  i  RO  (  3  )  ,  X<  500),  Y<  500  )  , 'J  (  50  0  )  ,  V  (  500  )  , 
1TH(1000) »STIF{ 1000*70) .  A  P  ( 1000) . EST I F( 6, 6 ) ,ECM( 3.3),EBM(3.6), 

_ 1 _ E_S  M(  3.6  J  .  M  Li _ NUMNP  « _ Mi  MEL  > _ NUMA.T_._KQ  D  El  5  Q  Q  )  «  NP(3 .1000)  . _ 

1  MAT (  1  000  ),  MBAND ,  NEQ .  M  *  L  M ( 6 )  .  I  MXM SH ,  I MES H{  50 0 ) ,  NCQRS(IOOO). 

1  IFSN( 20  ) ,NFP(20)  ,  N F SC , MFF , A P S (  1 00 0  )  , I H AR *  I NCOR ( 500  ) 

_ 1  .XKO.NCYC  ,MRP,DEPG.DEPE,NBP(20>  _ _ _ ...  _ 


C 

C  INITIALIZE  APPLIED  LOAD  VECTOR  AND  MASTER  STIFFNESS  MATRIX  AND  EC M 
£ 


DO  10  1=1  *  NE  Q 

AP (  I  )  =  0  •  0 

DO  10  J=1 , MBAND 

I  0 

ST  I F (  I  ,  J)=0.0 

DO  20  1=1.3 

DO  20  J=1 . 3 

20 

C 

c 

ECM (  I  .  J ) =0. 0 

FORM  ELEMENT  CONSTITUTIVE  MATRIX  (ECM)  IF  NUMA  T= 1  ) 

C 

IF(NUMAT.N£.l )  GO  TO  30 

EPR -PR ( 1 ) 

COM  =E ( I  ) / (  ( 1  . +EPR  )*  (1  .  —  2  •  * EPR  )  ) 
COM  1 =  CCM* (  1  .  — FPR ) 

CC  M2- CCM  +  EPR  _ _ 

ECM ( i , I )=C0M1 
ECM (2»2)=C0M1 


C 

tLM  l  1  *  S  J  -V.UI'Or' 

ECM ( 2, 1 )=COM2 

ECM ( 3 , 3 ) =E ( 1  )/(2.*(  1 • +EPR ) ) 

30 

C 

DO  45  M= l . NUMEL 

CALL  ELST IF( 1 ) 

C 

c 

ASSEMBLE  ELSTIF  INTO  MASTER  STIFFNESS  MATRIX 

DO  35  1=1.3 


12=2*1 


_ LM{  I  2 )  =  2*NP(  I  . MJ _ 

35  L M (  12-1  )=LM<  I  2 )  —  l 

DO  40  1=1.6 

I  I=LM{  I  )  _ _ 

DO  40  J  =  1  .6 
JJ  =  LM( J )- I  I  +  1 

J  F  ( JJ.LE.O)  GO  TO  40 _ 

ST IF( I  I  . J J) =  ST  IF(  I  I  .  JJ  )+EST  IF (  I  ,  J  ) 
CONTINUE 


40 
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c 

c 

AOO  GRAVITY  LOADS  INTO  AP  VECTOR 

00  45  1=2, 6, 2 
- 1 I=1_M.LL) _ 

45 

r 

API  I  I  )  =  A P (  I  I  )~WT 

CONT  INUE 

C 

C 

ADD  NODAL  LOADS  INTO  AP  VECTOR 

- DO  50  N= 1  , NU  MNP 

50 

N2=2*N 

API  N2 )=AP I N2 ) +V IN) 

API N2-1 ) =API N2-I ) +U I N ) 

* 

C 

c 

c 

MODIFY  STIFFNESS  AND  LOAD  VECTOR  FOR  DISPLACEMENT  BOUNDARY 

CONDI  tion; 

DO  100  N  =  1  , NU  MNP 

N2=2*N 

-  IF  (KOOE(N)-IO)  80,70,60 

6  0 

I  I  =  <N  2  —  l 

CALL  MODAFYI I I , N ) 

CAL'  M QELAfLYJ  N2.N  1 

70 

GO  TO  100 

I  I--N2-1 

CALL  MODAFYI I I , N ) 

80 

GO  TO  100 

I  FI KODEI N )  ,EQ  .0 )  GO  TO  100 

C  A  L  L  V,  0  D  A  FYI  N2.N) 

100 

C 

_ c _ 

CONTINUE 

_  OUTPUT  BODY  FORCES  FOR  LAST  CYCLE 

c 

IF  I  IHAR . NE. NCYC )  RETURN 

WR ITEI 6, 51 1 

51 

F ORM AT ( /// 24 F  TOTAL  NODAL  BODY  FORCES/ 

1  5X ,  5H  NODE , 1 0  X ,  8M  X  FORCE, 10X,  8H  V  FORCE/) 

_ DO  2005  NZ  =  1  »  NUMNP 

IY=2*NZ 

IX=  I  Y-l 

XFORCF=AD( IX) 

_ 52 

YFORCE=AP( l Y) 

WR I TE l 6 , 52 )  NZ , XFQRCE ,  YFCRCE 
_  FORM  AT  I 7  X ,  I3,10X,F1 0.1  ,8X,F10.  1) 

2005 

CONT INUE 

RETURN 

END 

MEMORY  REQUIREMENTS  000768  BYTES 
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SUBROUTINE  ELSTIF(KCP) 

C  THIS  SUBROUTINE  FORMS  THE  ELEMENT  STIFFNESS  MATRIX  (ESTIF)  OR 
C  ELEMENT  STRESS  MATRIX  (ESM)  FOR  THE  CONSTANT  STRAIN  TRIANGLE 

c : 

COMMON/ COM ING/  E ( 3  )  . FR ( 3 )  , RO ( 3 )  .  XI5C0).  Y(  500  )  .U ( 5 0 0 ) ♦ V ( 500 )  . 

1 TH  (  1000)  »STIF(  1000.70)  .AP(IOC-O)  »ESTIF(6»6)»ECM(3«3)»E  BMC  3.6). 

_ _ 1  ESM (3*6)  «  W  T  »  NU  MNP  <  NUVEL.  NUMAT,  KO0E(500)»  NP(3.1000)» 

1  MAT (  1 000 )  . MBAND .  NEU .  M,  LM(6>,  I M XM SK ,  I  MESH ( 50 0 )  »  NCORS(IOOO), 
i  IFSN( 20  )  ,NFP(20  )  ,  NF SC , MFP , APS (  1  00 0  )  ,  I H AR . 1 N COR ( 5 00  ) 

_ 1 «  XKO  »NCYC  «MBP .OEPG.OEPE. NB  F ( 2  0 ) _ 

C 

DIMENSION  AV(3).eV(3) 

DO  10  1=1.3 

DO  10  J  =  1  ,6 

_ LH _ FBMI  I.  J)=0.0 _ _ 

I =MP{ 1 , M) 

J=NP(2,M> 

_ _  K-NP  (  3  ,  M  )  .  .  .  _ _ _ _ _ _ 

C 

C  FORM  ELEMENT  DIMENSIONS 

c 

A V (  1  )  =  X(K  )  —  X ( J  ) 

A  V ( 2 ) =  X (  I )  —  X ( K ) 

AV  (  3  )  =  X  (  J  )  —  X  (  I  )  _ _ _ 

8 V ( 1 ) =Y (J)-Y(K) 

BV( 2 )  =  Y{ K )- Y(  I  ) 

_ BV( 3) =Y(  I ) - Y (  J  ) _ 

ARE  A2  =  AV ( 3) *BV ( 2 )-AV( 2  )*BV (3 ) 

IF  (TH(M).EQ.O.)  TH ( M ) = 1 . 

VOL  =  TH( M > *AR£A2/2  .  _ , _ : _ 

IF  (VOL.LE.O.)  GO  TO  75 
C  FORM  CONSTITUTIVE  MATRIX 

_ C _ 

IF  (NUMAT. EQ.l)  GC  TO  30 
NM-MAT (M ) 

EPR  =  PR  (  MO  _  .  ..  _ _ 

CO  M= E ( NM ) / (  (  1 . +EPR ) * (  1 .-2 • *EPR) ) 

COM  1 =COM*( 1  .-ERR ) 

_ CCLM2_S_C.  CaM  *_EP  R _ 

ECM ( 1  ,1  ) =  C  O  M 1 
ECM (2.2) =C0M1 

_ ECM  (  1  ,2  )=CCM2  _ _ 

ECM ( 2 . i ) =CCM2 

ECM ( 3,3)=E(NM)/(2.*(1.4EPR)) 

C _ 

C  FORM  ELEMENT  B  MATRIX  (EBM) 

C 

_ 30 _  DO  4  0  1  =  1.3 _ _ _ _ 

I  2  =  2* I 
I  2  M=2  *  I  —  1 

_ E  BM (  1 . I 2M)=BV( I 1/AREA2 _ 

EBM(2,I2)=AV( D/AREA2 
EBM( 3. I2M)=EBM(2« I? ) 

40  EBM(3. I2)=EBM( 1 , I2M)  _  _ 


. 


f 
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C 

C  FORM  ELEMENT  STRESS  MATRIX  (ESM) 

C 


-DQ S  0 


50 

C 

DO  50  J=1 ,6 

ESM(  I  ,  J)=0.0 

DO  5  0  K=l,3  _  _ _ _ 

ESM  (  I  ,  J  )  =ESM  (  I  ,  J  )  +  F.CM  (  I  .  K  )  *E8M  (  K  *  J  ) 

IFIKOP.EQ.?)  CO  TO  60 

c 

I  F ( NUMAT . EQ . 1  >  NM=l 

WT=VOL*RC( NM)/3. 

c 

c 

FORM  ELEMENT  STIFFNESS  MATRIX  (ESTIF) 

DO  70  1=1 .6 

OO  70  J=l ,6 

ESTIFI  I  .  J  >=0.0 

00  6  0  K=1  .3  .  _  _ 

60 

70 

ESTIFI I i J) =EST I F ( I . J ) +E E M ( K . I > *ESM ( K , J > 

EST I F { I , J )=EST IF ( I . J ) *VOL 

GO  TO  60 

75 

60 

V.RI  TEI  6, 1000  )  M 

CALL  EXIT 

RETURN 

C 

c 

1000 

FORMAT  (  1  HI  «  18H  VOLUME  OF  ELEMENT  ,14.  1 8H  IS  LESS  THAN  ZERO) 

c 

END 

.  MEMORY  REQUIREMENTS  000736  BYTES 
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SUBROUTINE  MOCAFY  (I,N) 

C 

C  THIS  SUBROUTINE  MODIFIES  KSTIF  AND  AP  IF  A  DISPLACEMENT  BOUNDARY  CCND 

_ C - IS  IMPOSED  IN  EQUATION  I  ASSOCIATED  WITH  NODAL  POINT  N _ 

C 

COMMGN/CCM INC/  E(3)  »PR(3)  »R0(3)  *  X(500).  Y (  50 0 )  ,U ( 5 0 0 )  , V ( 500 )  . 

_ 1 TH (  1000)  ,ST IF( 10  00  *70 ) . AP(  1 00  0  >  • EST I F (6 ,6 ) »ECM<  3  *3 ) ,E8M( 3*6) • 

1  ESM(3»6),WT,  NUVNP,  NUMEL.  NUMAT,  K.0D£(500)  .  NPI3.1000)* 

1  MAT(  1000  )  . MBAND ,  NEO,  M,  LM ( 6 )  ,  I M X M SH ,  I MES H ( 50 C ) .  NCCRS(IOOO)* 

_ 1  I  ESN  I  20  )  .  NFP  (20  )  .  NF  SC  .  V.  F  P  .  A  P  S  (  1  0  Q  0  )  .  I  HAP  .  I  NCOR  (  5QQ  i _ 

1  *  XKG  »  NC Y  C»MBP»DEPG*DEPE»NHP(20) 

C 

_ _ D  I  SP  =  U  (  N  )  .  _ 

IF{ ( I-2*N) .EQ.O )  D I SP= V ( N ) 

C 

_ DO  50  J=2. MBAND _ 

I  L=  I  +-  J  —  1 
IU=I-j+l 

_  IF ( IU.LE.OIGO  TO  10  _ 

API  I U)-APl  IU )  -  STIFI  IUt J )*DI SP 
ST IF{ IU  . J)=0.0 

_ LO _ IF(  IL.GT.NEQ)  GO  TO  50 _ 

APilL)=  API  1 L ) —  STIFI  I  • J ) *D I SP 
STIFI  1  ,  J)=0*0 

50  CONTINUE  _  _ _ _ _ _ _ 

API  I ) -D  ISP 
STIFI I » 1 )=1 . 0 

_ RETURN _ 

END 

..MEMORY  REQUIREMENTS  0002DE  BYTES  - 
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C 

C 

c 


SU8R0UT I NE  RANDl 

This  subroutine  solves  for  displacements  USING  a  GAUSSIAN  elimination 

-TECHNIQUE  FOR  < 


CCMMCN/CCM INC/  E ( 3 )  , PR (  3 ) . RO { 3  )  ,  X(500).  Y ( 500 )  . U ( 50 0 ) . V l 500 ) , 

- 1THI1000),  A(1000,70).  B ( 1 00 0 ) . EST I F ( 6 , 6 ) , ECM( 3 . 3 ) ,E BM ( 3 , 6 > , 

1  ESM(3,6),WT,  NU*NP,  NUVEL.  NUMAT,  KODE(500),  NP( 3,1000) , 

l  M AT ( 1 000 ) *  MM  .  NN  •  M«  LM(6),  I M XM SH ,  I  MESH { 50 0 ) •  NCORS(IOOO). 

- - 1 - 1  ElSNJ  P  (  20  )  .  NFSC  .  MFQ  .  APS  (  1  000  )  .  I  HAR  .  I  NCOR  <  5  UP  ) _ _ 

1  iXKOiNCYC »MBP  »DEPG,DEPE, N8P(20) 

C  TRI ANGULARI Z£  AND  REDUCE  RIGHT  HAND  SIDE 
_  NL=NN-MM+1 

NM=NN  — l  "  . '  " 

MR  =  MM 
X 

DO  100  N=1,NM 

IF  ( A(N» 1 ) .LE.O. )  GO  TO  700 

_ BN=B(N) 

BIN) =BN/A { N ♦ 1 ) 

IF  (N.GT.NL)  MR=NN-N+1 
X 

DO  100  L=2 • MR 

IF  ( A(N.L) .EC.O. )GO  TO  100 

_ C=  A(N,L)/A<  N. I  ) 

I-  N+L-l  - 

J-  0 

_ DO  50  K  - L  .MR  _ . _ 

J  =  J  +  1 

50  All  «  J ) =  A(  I  ,  J  )-C*AC N,K) 

_ 3  (  I  )  =  B(  I  )  —  C  *  8  N  _ _ 

A(N,L) =C 
100  CONTINUE 

X _ _ _ 

C  BACK  SUBSTITUTE 
C 

_ I  — NN  _  .  _  _ 

B(NN)=0( NNl/AINN, 1 ) 

DO  500  N-l.NM 

_ 1  =  1-1 _ _ _ _ 

I F  (N.LT.MM)  MR  =  N+l 

DO  600  J  =  2 1  MR 

_ K-UJ-i  _ 

600  0(1)  =B( I )—  A(I,J)*B(K)  * 

IF(  IHAR.NE.NCYC  )  RETURN 

X _ _ _ _ 

C  OUTPUT  DISPLACEMENTS  FCR  LAST  CYCLE 

C 

WRI T  E ( 6 ♦ 209 1  )  1  FAR _ 

WRITE! 6,2001  ) 

2091  FORMAT! 1H1.  11H  CYCLE  NO.=,  13/) 

2001  FORMAT! _ «  NODAL  DISPLACEMENTS  »//» _ NODE  X-  D  I  SPL  ACE 

1MENT  Y-DISPLACEMENT*//  ) 

DO  2005  NZ=1 1 NUMNP 
I  Y  =  2*NZ 


' 


4 
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I X= I Y- 1 
XDISP=8< IX) 

YOISP=E( IY) 

J2  0Q-5 _ WRITFt^^nnft . 

2006  FORMAT  (  I  9 , E2 1 . 5 ,E 2 0 . 5 ) 
RETURN 

.700  _ WRITE  (6  .2000)  N  _ 

CALL  EXIT 

2000  FORM  AT (  1  HO .  81H  Z! 
LNGULAR.IZED  MATR! 


:rq  cr  negative  element 

.  15  ) 


ON  MAIN  DIAGONAL  CF  TRIA 


END 
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SUBROUTINE  STRESS 
C 

C  THIS  SUBROUTINE  FORMS  THE  ELEMENT  STRESS  MATRIX  I E SM ) . M UL T I  PL  I ES  BY 

- C - 1  HE — FLF  MEN  I — DISPLACE  WENT — VECT  OR.  XELDLEPJ  AND  _REjCQEi2_S_THE__ST  RESSFS  IN 

SIGEL.  IT  THEN  COMPUTES  THE  PRINCIPAL  STRESSES  AND  D I RECT I C NS ( S I GP ) 

- CCMMGN/COMING/  E ( 3 ) « PR ( 3 ) .R0{3 ) .  XI500),  Y { 5 00 ) . U ( 5 0 0 ) . V ( 500 ) . 

t TH{  1000  )  ,ST  IFl  1 000 ,70)  , API  1000)  .ESTI F( 6,6)  ,ECMl 3,3) .EBMl 3.6J  , 

1  ESM(3»6),WT,  NU  MNP.  NUMEL.  NUMAf  »  KUDEI500  ).  NP ( 3 , 1 0  0  O )  » 

- - - 1 — WAILLQ.O  Q_L«  MfLANQ  • — NtLQji — Mj _ LMJ_6Jj _ I MXM5H ,  I ME  SHI  500 ) ,  NCGRSIIOOO). 

1  IFSN I  20  )  , NFP  l  20  )  ,  NFSC, MFP, APS (  1000  )»  I H AR ,  INCOR ( 500  ) 

1  .XKO.NCYC  »MBP,OEPGf OE PE  .NBP120) 

- DIMENSION  SIGELI  1000,3)  ,  .....  _  ELD  I  SP  16),  KOUNT  (  500  ) 

COMMON/ETS/SIGPI 1000.4) 

c 


_ DO  5  N=  1  «  NUMNP _ _ _ _ 

KOUNT (N)=0 
5  CONTINUE 

_ C _ _  .  _ _ _ 

DO  100  M  =  1  ,  NUMEL 
C 

- C - COMPUTE  ELEMENT  DISPLACEMENTS _ 

C 

DO  10  1=1,3 

_ 12  =  2*1 

LMI2  =2*NPII,M)  '  "  '  1 

ELD  I SP (  12) = AP ( LM I  2 ) 

- U2 _ ELDISPI  12-1  )  =  AP  I  LM  12—1  ) _ 

C  ~ 

C  COMPUTE  ELEMENT  STRESSES 
_ C 

CALL  ELSTIFI2) 

C 

_ DO  20  1  =  1,3 _ _ 

SIGELI M, I ) -0 .0  ~~ 

DO  20  J=1 ,6 

_ 2.0. _  SIGELI M,  I  )  =S  I  GEL  I  M  ,  I  )  *ESM  I  I  ,  J  )  *ELDI  SPI  J) 

C  FIND  AVERAGE  X  AND  Y  OF  ELEMENT  M 

J1 =NP I l , M ) 

_ J2.  =  N  P  I  2  ,  M  ) _ 

J3=NP13.M) 

X A VE=0 .3333*(XIJ1)+XIJ2)+X(J3)) 

_ YAVE=0. 3333*1  Y{  J1  )  *Y  I  J2)+YU3)  )  _  _ 

C 

C  CALCULATE  THE  AVERAGE  PRINCIPLE  STRESS 

C 

APS  I M )=0 .3  13  3* l-ROI 3 ) *1 CEPG-YAV^  )  *1  1 . +2.4XKO  )*(  SIGELI v,  i  )  + 
1  SIGELI «  ,2 )  ) *  I  1  .+PRI3  )  )  ) 

_  .  _  S  I  GELt M,  1  ) =S  I  GEL  I m ,  l  )  - t  t DLPG-YAVE ) *00  I  3 ) *XKU ) 

SIGELI M, 2  )=  S I GELI M , 2  )- I  I UEPG-Y AVE  )  *RO I  3)  ) 

B  3  Cf NT  I NUE 


C  COMPUTE  ELEMENT  PRINCIPAL  STRESSES  AND  DIRECTIONS 
C 

SIGM=  (S  I  C-EL  I  M,  1  )4SIGEL(M,  2)  )  /2. 


prr- 
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S  I  GD2=  (  S  l  Gt  l.(  M,  l  )  -S  IGEL  (  M  .  2  )  )  /2  . 
RAO  =SORT ( S I GD2**2  +S  IGEL ( M , 3 ) **2 ) 
S  IGP ( M , 1 ) =S 1 GM+PAO 
-S.IGPlK.21-_ 


-S1GNL-RAD, 

RAO 


_1  00 
c 

r 


S IGP ( M , 3  )  = 

CONT I NUE  J  =  O,S*57’29570*ATAN2 (SIGEL(M»3>  .SIG02) 
IF ( IHAR.NE.KCYC)  RETURN 


C 

c 


— EIND - MAXJMLLM _ ELEMELN  T  .S.TRESS£_S__ 

OUTPUT  STRESSES  FOR  LAST  CYCLE 

_ SIG1=0.0 

SIG2=0.0 
SIG3=0 .0 
Ml  =  0 


M2  =  0 
M3=0 


_ _ OO  130  M=1,NUM£L 

IF  ( S I GP { K  .  1  )  ,LT . SI Gl  )  GO  TO  115 
SIG1=SIGP{ M, 1 ) 

M1  =  M 

115  IF  ( S IGP ( M  ,2 )  .GT .S IG2 )  GC  TO  120 
SIG2=SIGP(M.2 1 
_  M2=M 

120  IF  { SIGP( M,3)  .LT.SIG3)  GC  TO  130 

SIG3  =  S  I  GP ( M.3  ) 

M3-  M 

130  CONTINUE 

WR l TE(6 ,2019  )  IMAP 

_ tePIT £(6,2020  )  _  _  _ 

WRI  TE  (6 .2030  )  (  M  .  (  S  I  GP  (  V  .  I  )  .1  =  1  .4)  ,APS(M)«M=1,  NUMEJ-  ) 

WRI TE (6 . 20  40  1  (3IG1  .Ml  .SIG2.M2.SIG3.M3) 

JL3JL1 _ CONTINUE _ _ 

C  ~~  ' 

C  FORMAT  STATEMENTS 

C _ _ _ _  _  _ _ _  _ 

2019  FORM  AT (  1H 1  .  11H  CYCLE  NG.=  ,  13) 

1306  FORMAT!  15. 3F1 0.3  ) 

-1.3  OQ  FORM  ATII4.F13.2.F12.2.F12.3) _ 

2000  FORMAT ( 1H 1 , l OX ,  21H  X-Y  ELEMENT  STRESSES  ///  40H  ELEM  SIGMA  X 

1  SIGMA  Y  SIGMA  XY  //) 

2010  FORMAT! 14, FI  1 .3.2F12.3)  _ _ 

2020  FORMAT!  /10X.27H  ELEMENT  PRINCIPAL  STRESSES  /// 

1  55  H  ELEM  SIGMA  I  SIGMA  II  MAX  TAU  SIGMA  I  DEG  , 

_ 1  9  X  .  5H  APES  //) _ _ 

2030  FORMAT (14* F 11.3* 2F12 .3*  F14.3.F16.1) 

2040  FORMAT ( 1  HI  . 

1  27H  MAXIMUM  PRINCIPAL  STRESS  =  .F11.1.I9H  AND  OCCURS  IN  ELEM, 16// 

2  27  H  MINIMUM  PRINCIPAL  STRESS  =  ,Fl  1.1.19H  AND  OCCURS  IN  ELEM, 16// 

3  27H  MAXIMUM  SHEAR  STRESS  =  .F11.1.19H  AND  OCCURS  IN  ELEM, 16) 

_ ENJ2 _ _ _ _ 
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C 

C  FLOT  PROGRAM 

C 

_ _ _ _ 

C  THIS  A  PROGRAM  WH  I  CH  TAKES  PUNCHED  OUTPUT  FPCM  PROGRAM  FPM5 

C  WHICH  CONSISTS  OF  ELEMENT  AND  NODAL  INFORMATION  AND  THE 

c _ POTENTIAL  AT  THE  NODES  AND  USES  I  T_  TO  OUTPUT  AND  PLOT  THE  DESIRED 

C  ECUIPOTENT I ALS 

C  THESE  EQU  IPOTENT I ALS  CONSIST  OF  STRAIGHT  LINE  SEGMENTS  BETWEEN 

_C _ THE  INTERPOLATED  POINTS  _ 

C 

c  •  the  FIRST  CARD  MUST  CONTAIN  INFORMATION  TO  IDENTIFY  THE  PLOT 

C  _  AND  WILL  E-:  THE  UNLY  IDENTIFYING  INFORMATION  WRITTEN  ON  THE  PLOT 

C 

C  CARD  2  MUST  CONTAIN  THE  FOLLOWING  INFORMATION:  THE  NUMBER  OF 

_C _ ECU  I  POTENT  I ALS  TO  bE  FOUND,  THE  NUMBER  OF  NODAL  POINTS.  THE _ 

C  NUMBER  OF  ELEMENTS,  THE  MAXIMUM  Y  VALUE  HOP  THE  F.E.  GRID 

C  (THE  PROGRAM  CONVERTS  THE  GRID  CO-ORDINATFS  INTO  PAGE 

c _  CC-CRDINATES - THE  GRID  SHOULD  BE  NEAR  THE  AXIS  IN  +X-Y  QUADRANT) 

C  AND  THE  DESIRED  HEIGHT  OF  THE  PLCT(THIS  SHOULD  3E  LESS 

C  THAN  25  INCHES) - FC RM A T ( 3  I  5 , 2F l 0 . 2 ) 

_C _ 

C  CARD  3  MUST  CONTAIN  THE  EQU IPOTENT I ALS  CESIREO 

C  - FORMAT  (FI0.2) 

C _ _  _ _ _ _ _ _  _ 

C  THE  CARDS  WHICH  FOLLOW  MUST  CONTAIN  INFORMATION  FROM 

C  PROGRAM  FRVE-  ONE  CARD  FOR  EACH  NODE  DESCRIBING  ITS  POSITION 

SL _ AND  POTENTIAL  FOLLOWED  3Y  THE  ELEMENT  I NF OR m A T I  ON  CARDS _ 

C 

C  FOR  EVERY  EQU I PO T ENT  I A L  THE  PROGRAM  TAKFS  EACH  ELEMENT  IN  TURN 

C _ AND  CHECKS  FOR  THE  DJES  IRED__EQU  IPOTENT  IAL _ 

C  IT  THEN  INTERPOLATES  AND  STORES  THESE  VALUES 

C  'a  MAXIMUM  Y  VALUE  IS  FOUND  AFTER  WHICH  THE  CLOSEST  POINT  TO 

c _ this  maximum  y  is  obtained  and  then  the  closest  point  to  the _ 

c  PRECEEDING  POINT  AND  SO  ON 

C  NOTE -  THIS  PROCEDURE  WILL  NOT  WORK  IN  ALL  CASES 

C _ _ _ _ _ 

C  THE  PROGRAM  HAS  DIMENSION  STATEMENTS  FOP  25  CONTOURS, 

C  5 C 0  NODES,  900  ELEMENTS  AND  99  INTERPOLATED  VALUES  FOR  EACH 

C _ CONTOUR _ _ _ 

C 

C  I  MENS  I GN  EPC(25),PSI(500)»XP(99»25)  ,YP(99«25)  »  XP  M  (  99,25)  , 


_ 1  YPM ( 99 .25 ) ,  X <5 00  ).  Y(500  ),  R ( 50 0  ) , DA T A (  40 0 0 ) 

DIMENSION  I  A ( 9  9 )  ,NP(4,900)  ,ICHK(50  0.2) 

INTEGER  HED (10) 

_ CALL  PLOTS ( C AT  A (  1  )  ,  1  6000  ) _ 

C  READ  INPUT 

RE  AC ( 5 , 4  )  HED 

_ RE  AD ( 5 , 9 )  NEPC  ,  NUMNP . NE , YMAX , DH _ 

READ(5,e)  ( EPC  (  I  )  ,  1  =  1  ,NF.PC) 

READ(5»7)  (  I  ,X(  I)  ,Y(  I  ),R(  I  )  , PS  I (  I  ),  1  =  1, NUMNP) 
_ DO  70  1  =  1  .NF _ 

7  FORMAT ( 3X ,  I 3»T  55,2F3.2. T 4 1 , 2F  7.2) 

6  FORMAT ( 12X . 51 3 ) 

_ REAP( 5 ,6)1  , i  ( N P(  J .  I  )  . J= 1 , 4) _ 


■ 


C3 


'»  COMP  ILfcfi  MAIN 


03-02-72  10:32.40  PAGE  0002 


70  CONTINUE 

5  FCPMAT ( 1H1 , 1 RA4 ) 

Jt_ F.CHMAT  (  GF  10.  I  ) 

9  FORMAT  (3  IE, 2F10.2) 

4  FORMAT (  18A4  ) 

_  *RITE( 6.5)  FED _ 

D032  L  =  1  »  NEPC 
I  A  (  L  )  =  1 
MCT=Q 

DC  10  1  =  1  , N  E 

DC  10  J=1  ,4 

_ K  =  NP<  J  ,  I  )  _ 

JP=J+I 

IF (NP ( 3 , I  ) . EC ,NP< 4 , I  )  )  GO  TO  81 

CO  TO  80 

81  I F ( J . EG • 3  )  GO  TO  10 

80  I F ( J • EG • 4  )  JP  =  1 

_ KK-NP  (  JP  ,  I  )  _ _ _ 

A=PS I ( K )-EPC(L ) 

P=PSI (KK)-EFC (L) 

IF((A*B)«C-T  .  (  C«0)  )  GO  TO  10 

C  CHECK  TO  SEE  IF  THE  INTERPOLATION  HAS  ALREADY  BEEN  DONE 

C  IN  AN  ADJACENT  ELEMENT 

MCT=MCT+l  _ 

ICHKI MCT , 1 ) =K 
I  CHK ( MCT ,2) =KK 
T  C  AT  =  0 

DO  77  JX= l , MCT 

IF (  ICHK ( J X ,  1  )  . NE • KK )  GO  TO  77 

_ IF ( I C HK ( J X. 2 ) . tO  *K) _ I  CAT  =  1 _ 

77  CONTINUE 

IF( ICAT.EO. I )  GO  TO  10 
I  1= I  a ( L  ) 

C  INTERPOLATE 

XP(  II  »  L  )  =  X  l  K )  +  ( (  (  X(KK)  -X{K))*(PSI(K) -EPC  ( L  )  )  )/  { PS  I  (K ) -PS  I  < KK  )  >  ) 

_ _  YP  (  I I  ,  L  )  = Y { K  )  +  ((( Y ( KK )  -Y ( K  )  ) »(PSI ( K) -EPC ( L  )  )  ) / ( PS  I  ( K ) -PS  I ( KK )  )  ) 

I  A  (  L  )  =  I  A  (  L  )  +  l 

10  CONTINUE 

IF (  I  A (L ) .LE.2 )  GO  TO  39 

GO  TO  38 

39  WRITE(6»1C1  )  EPC(L) 

_ GO  TO  32 _ 

38  CONTINUE 

101  FORMAT (• 0 56H0NE  OR  ZERO  INTERPOLATED  VALUES  FOUND  FOR  EOUIPOTEN 

_ 1  T I AL= .  F  6 . 1  ) _ 

C  FIND  THE  MAXIMUM  Y  OF  THE  E CU I PO TEN T I  A L  VALUE  IN  QUESTION 

YP  M ( 1  , L  )=-!  0000 . 

_ I  I  =  I A  (  L  )—  1 _ _ 

DO  20  1=1,11 

IF(YP{ I  ,L  )  .LT.YPMl 1 ,L)  )GO  TO  20 

_ YP  M (  1  ,L ) =YP  (  I  ,  L  ) _ 

XPM( 1,L)=XP(I,L) 

IC=I 

20  CONTINUE  _ _ _ 


* 
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C 

11=11-1 

IF ( IC .GT . I I )  GO  TO  40 

MOVE  ALL  VALUES  FROM  THE  MAXIMUM  FOUND  TO  THE  LEFT  ONE 

DO  40  I  =  IC .  T  I 

LOCATI ON 

40 

YP (  I ,L ) = YP (  1  +  1,0 

XP< I ,L)=XP( 1  +  1,0 

CCNT I NUE 

C 

FIND  NEXT  CLOSEST  POINT 

DC  31  I M=  1  »  I  I 

J=I  I-I M+ 1 

EC={(YPM(IM  ,L  )  — YP(  1 .L )  ) **2+ ( XPM(  I M  . L )- XP <  1  .  L ))**?>** 0 . 5 

IC=1 

IFIJ.LT.  2)  GC  TO  30 

DC  30  1=2, J 

0O=((YPM(IM  ,L)-YP(  I  ,L)  )**2+(XPM{  [M  , L ) - XP (  I  , L  )  ) ** 2 ) * * 0 . 5 

IF(BO.LT.OC)  GO  TO  28 

28 

GC  TO  30 

IC=I 

BC=6D 

30 

CONTINUE 

YFV(  I M+l  »L)=YP(  I C  »  L ) 

XPM( I M+ 1 ,L ) =XP ( ICiL) 

C 

I  J=I  I-I M 

I F ( IC.GT. IJ )  GO  TO  50 

MOVE  ALL  VALUES  TO  THE  LEFT  ONE  LOCATION 

CO  50  I  =  IC.  IJ 

YP (  I »  L ) =  YP<  1  +  l  ♦  L  ) 

XP(  I.L)=XP(  T  *  1  •  L  ) 

50 

3  l 

32 

CONTINUE 

CONTINUE 

CONT  IN'JE 

DC  15  1=1 .NEPC 

I  I  =  I  A (  I  )~  1 

IF(  I  I  .LE.  1  )  GO  TO  15 

• 

WRITE(6,6fc)  EPC ( I ) 

¥RITE( 6,65) 

WRITE<6,67)  ( XPM ( J  »  I  >  ,  YPM ( J  ,  I  )  , J= 1  ,  I  I  ) 

15 

66 

65 

CONT I NUE 

FORMAT (’O'*  32H  CO-ORDINATES  FOR  E OU I POTEN T I AL= . F 1 0 . 2 ) 

FORMAT ( 15X, 2F  X.15X.2H  Y) 

67 

FCRMAT(9X,F10.2,7X,F10.2) 

YP0S=DH+3  .0 

CALL  SYMBOL { 1 .0 , YPOS . 0. 20, HEO. 0.0.72) 

C 

PLOT  GRID 

SF=YM AX/DF 

DC  23  1=1 , NUMNP 

23 

X (  I  )=Xl  I  )/SF+1.0 

Y  (  I  )  =  Y  (  I  )/SF  41  .0 

CCNT I NUE 

DO  24  1=1 »NEPC 

I  I  =  I  A (  I  )-  1 

I F ( I I .LE . 1 ) GO  TO  24 

OO  24  J=  1  ,  I  I 

XP  M ( J ,  I  )=XPM( J,  I  )/SF ♦ 1  .0 

YPM  J  ,  I  )  =  YP  M  (  J  »  1  )/SF  *1  .0 

. 
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24 

C  C  NT  I NUE 

CO  11  J=1.NE 

K  =  NP< 1  ,  J) 

CAll  PI  flT  (X  f  K  )  •  Y  (  K  )  .  3  J 

DO  1010  1=2,5 

L=I 

IF ( NP(  3 , J )  .  EQ. NP(  4  ,  J  )  )  GO  TO  91 

91 

PI  l 

GO  TO  811 

IF (  I . EO  .  4  )  GO  TO  10  10 

TFfI.FO.5)  L  ~ 1 

10  10 

1 1 

K=NP(L , J  ) 

CALL  PLOTIX (K) ,Y{K) ,2> 

CONTINUE  -  .  -  -  _  —  _ 

CC  NT INUE 

DC  16  1  =  1  , N  EPC 

T  I  =  T  A  (  I  )  -  1 _ — - 

IF (I  I  .Lt. 1  )  GO  TO  16 

CALL  PL OT (  X  P  V  (  1  *  I  )  *YPM{ 1  «  I  ) t3) 

rnirm.i^i.TT  ....  _  - 

100 

1  9 

CALL  PLOT  (  XPM  (  J  ,  I  )  ,  YPM  (  J  »  I  )  .  2  ) 

CONTINUE 

PX-XPM  (  I  I  ,  I  )— 0  •  l  5 

1  6 

PY  =  YPM(  I  I  ,  I  )-0 .30 

CALL  NUMBER (PX.PY.O. 15,EPC( l ) ,0.0,1 ) 

CONTINUE  -  - 

CALL  PLOT (0 .0 ,0. 0 .999 ) 

END 

]RY  REQUIREMENTS  015340  BYTES 
l  RC  =  0 
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ITEM- 3  CONTINUED  ITEM- 3 

UPPER  BARREL 
(Stainless  Steel) 
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ITEM  -  4 

LOWER  RUBBER  PACKING  RETAINER 
(Stainless  Steel) 


0.125 
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ITEM  -  5 

UPPER  RUBBER  PACKING  RETAINER 
( Stainless  Steel ) 


0.125 


035 


CUT  THROUGH 


1.25 


o 

Cm 
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ITEM-6 

UPPER  CAP 
(Stainless  Steel) 


1.45 


h 


To  fit  1.25 


I.D.  Machined 
so  as  to  fit 
over  O-Rings 
&  Rubber  mem- 


o 

o 


brane. 

ITEM  -  7 

SEALING  MEMBRANE  RETAINER 
( Stainless  Steel ) 
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ITEM  -  8 


ITEM-  10 

PORUS  TIP'RETAINER 
(Stainless  Steel) 
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